Accepted by the ApJ 



Regulation of Star Formation Rates in Multiphase Galactic Disks: Numerical 
Tests of the Thermal/Dynamical Equilibrium Model 

^ Chang-Goo Kim^'^, Woong-Tae Kim^'^'^, and Eve C. Ostriker^ 

^ Center for the Exploration of the Origin of the Universe (CEOU), AstrononiAj Program, 
Department of Physics & Astronomy, Seoul National University, Seoul 151-742, Republic of Korea 

bJO 

^ ^Department of Physics & Astronomy, FPRD, Seoul National University, Seoul 151-742, Republic 

^ of Korea 

^ ^Institute for Advanced Study, Einstein Drive, Princeton, NJ 08540, USA 

^ "^Department of Astronomy, University of Maryland, College Park, MD 20742, USA 

o 

kimcgOastro . snu . ac . kr , wkimOastro . snu . ac . kr , ostriker@astro.umd.edu 
O ABSTRACT 

I ^1 We use vertically-resolved numerical hydrodynamic simulations to study star forma- 

tion and the interstellar medium (ISM) in galactic disks. We focus on outer disk regions 
^ where diffuse H I dominates, with gas surface densities S = 3 — 20 Mq pc~^ and star- 

00 plus-dark matter volume densities psd = 0.003 — 0.5 M© pc~^. Star formation occurs in 

very dense, self-gravitating clouds that form by mergers of smaller cold cloudlets. Tur- 
bulence, driven by momentum feedback from supernova events, destroys bound clouds 
^ and puffs up the disk vertically. Time-dependent radiative heating (FUV from recent 

^ star formation) offsets gas cooling. We use our simulations to test a new theory for 

^ self-regulated star formation. Consistent with this theory, the disks evolve to a state 

K>" of vertical dynamical equilibrium and thermal equilibrium with both warm and cold 

phases. The range of star formation surface densities and midplane thermal pressures 
^ is SsFR ~ 10""^ - 10-2 Mq kpc-2 yr"^ and Pth/^B ^ 10^ - 10"^ cm'^ K. In agreement 

with observations, turbulent velocity dispersions are ~ 7 km and the ratio of the 
total (effective) to thermal pressure is Ptot/-Pth ~ 4 — 5, across this whole range (pro- 
vided shielding is similar to the Solar neighborhood). We show that Ssfr is not well 
correlated with E alone, but rather with Sy^/Ogd) because the vertical gravity from stars 
and dark matter dominates in outer disks. We also find that Ssfr has a strong, nearly 
linear correlation with Ptot; which itself is within ~ 13% of the dynamical-equilibrium 
estimate -Ptot.DE- The quantitative relationships we find between Ssfr and the turbu- 
lent and thermal pressures show that star formation is highly efficient for energy and 
momentum production, in contrast to the low efficiency of mass consumption. Star 
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formation rates adjust until the ISM's energy and momentum losses are replenished by 
feedback within a dynamical time. 

Subject headings: galaxies: ISM — galaxies: kinematics and dynamics — galaxies: star 
formation — method: numerical — turbulence 



1. Introduction 



Large-scale star formation rates in galaxies are observed to correlate with both the gaseous 



and stellar content, and with the galaxy's gravitational potential well (e.g. 


Ryder &; Dopita 


1994 


Kennicutt||1998 


Wong & Blitz||2002; Boissier et al.||2003; Salim et al.||2007; Leroy et al.||2008 Bigiel 


et al. 


2008 


2010 


2011 


Genzel et al. 2010 


Daddi et al. 2010 


Shi et al. 


2011 


). Empirical fits in 



disks often adopt power- law ("Kennicutt-Schmidt") forms for the relationship among the surface 
density of star formation Ssfr, the surface density of gas S, the surface density of the old stellar 
disk Eg, and the orbital angular velocity fi. 

From the "supply side" point of view, gas represents the fuel for star formation, and the stellar 
disk and dark matter halo help to define dynamical timescales within the interstellar medium (ISM) 
that could affect how rapidly gas collects and collapses: the galactic orbital time, the vertical 
oscillation period and flow crossing time, and the gravitational free-fall time. Power laws naturally 
arise if the star formation rate is proportional to the ratio of the gas content and one of these 
dynamical times. The observed timescale for gas to be converted to stars, tsF,gas = 5]/Ssfr is, 
however, generally quite long compared to these dynamical times. Together, the empirical results 
present a picture of star formation that is sensitive to both fuel supply and ambient environmental 
conditions, and that has low apparent efficiency. 



In recent work, Ostriker et al. (2010) (hereafter OMLIO) and Ostriker & Shetty (2011) (here- 
after OSll) have argued that star formation rates respond to demand, as well as supply. Main- 
taining an equilibrium state in the ISM requires constant inputs of energy and momentum, and 
contributions from star formation are critical. Star formation can be self-regulated via feedback, 
in such a way that supply and demand match within the ISM: heating balances cooling, pressure 
balances gravity, and turbulent driving balances dissipation. The theory of OMLIO and OSll 
proposes that observed star formation rates can be understood as a response to the needs of the 
ISM. Because each massive star injects so much energy, only a relatively modest star formation 
rate (implying a long tsF.gas) is necessary. From the point of view of energy and momentum sources 
and sinks, star formation is in fact quite efficient. 

To see why feedback is vital, it is key to consider the internal thermal and dynamical state 
of the ISM, rather than just integrated properties. The internal vertical dynamical time tdyn oc 
{Gptot)~^^'^ , for Plot the total (gas + stellar) density, depends on the thicknesses of the gaseous and 
stellar disks. In particular, the contribution from gas gravity alone gives Y^/tdyn oc Y?!"^ j H^l"^ . The 
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gas disk thickness H depends (linearly or quadratically) on the vertical velocity dispersion of the 
gas, which includes both thermal and turbulent termsj^ Because thermal energy is radiated away, 
and turbulent energy is dissipated (in shocks and shear layers) on timescales ^ tdyn ^ *SF,gas) the 
internal energy must be continuously replenished. Young, high-mass stars restore this energy and 
preserve the life of the ISM. If star formation feedback were entirely absent and the only heating 
source were the cosmic background radiation, tjyn would drop by nearly two orders of magnitude, 
with a corresponding (or greater) increase in SgpR- 

OMLIO and OSll, considering respectively mid-to-outer disks and central starburst regions, 
showed that observed star formation rates are quantitatively consistent with analytic predictions 
that follow from imposing thermal and dynamical equilibrium in the diffuse ISM. OSll also pre- 
sented initial results of numerical simulations that include turbulent driving associated with star 
formation, confirming the analytic theory for molecule-dominated regions. Additional results from 
simulations in the starburst regime will be presented in Shetty & Ostriker (2011, in preparation). 

In this paper, we use time-dependent numerical simulations to test the OMLIO theory (and 
extensions based on OSll), for the outer-disk regime where the ISM is dominated by diffuse atomic 
gas. A crucial aspect of our simulations is that we vertically resolve the disk (our grid scale is 1 pc). 
We shall show that, as assumed by OMLIO, thermal and vertical dynamical equilibrium are both 
satisfied in our numerical models. We shall also show that feedback from star formation is largely 
responsible for sustaining both the thermal and turbulent pressure (and energy) in the atomic ISM. 
We numerically calibrate the yield relation between Ssfr and the thermal and turbulent pressures in 
the diffuse ISM, demonstrating that near-linear relations hold for both Pth and -Pturb- By combining 
these feedback relations with dynamical equilibrium, we show that Ssfr depends nearly linearly 
on the weight of the diffuse ISM (i.e. the dynamical-equilibrium pressure -Ptot.DE ~ -fth + ^turb)- 
The correlation between SgpR and Ptot (or Ptot,DE) is stronger and more general than other star 
formation relations that are commonly cited. 

In addition to testing the thermal/dynamical equilibrium theory of star formation, our nu- 
merical models allow us to address a number of interesting issues related to observations of diffuse 



atomic gas in the Milky Way and external galaxies (Dickey et al. 1990: Braun 1997; van Zee &; 



Bryant|[T999t [Heiles Trda^Idl[2003l [Young et all[2003l [Petric Ru"pml[2007l [Dickey et aL][2009 



Kalberla & Kerp 2009). These observations show that (1) turbulent velocity dispersions are typ- 



ically ^ 10 km s~^, relatively independent of location or star formation rate; (2) both cold and 
warm atomic gas are pervasive, in proportions that appear relatively independent of location; (3) 
the thermal pressure is a small fraction of the total pressure. Our numerical results are consistent 
with these observations, and can be understood based on the thermal/dynamical equilibrium model 
with energy and momentum feedback from star formation. 



^ In this work we neglect the magnetic term, which is hkely to be small (see below) but would provide a minimum 
vertical support in the limit of vanishing turbulent terms. 
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Our numerical models are idealized in that they represent a local patch of unmagnetized gas 
in a featureless disk where star formation is primarily responsible for the injection of thermal and 
kinetic energies. Thus, in this paper we do not capture the potential consequences of galactic 
structural features and certain instabilities that may affect ISM dynamics and star formation. 
The ISM surface density averaged over ~ kpc scales can be significantly affected by large-scale 



gravitational instability (e.g. 


, Wada & Norman||1999 


2007 


Kim & Ostriker ||2001 


2002 


2007 


Li et 


al. 2005; Tasker & Bryan 20 


D6 Tasker & Tan|2009 


Tasker 


2011t|Bournaud et al.||2007 


Bournaud 


& Elmegreen 


2009; Hopkins et al.||2011), spiral arm compression (e.g. Kim & Ostriker 


2002 


2006 


Shetty & Ostriker|2006; Kim et al.|2008 2010| Dobbs Hi Bonnell|2006 , 2008; Dobbs et al.|2008 


2011 


Wada & Kod; 


I 2004 


Wada 5 


!008 Wada et al. 2011) 


and Parker instability (e.g. 


Basu et al. 


1997 


Kim et al. 1998 2001 


L 2002 


Mouschovias et al. 200! 


) ) . Since the timescales to collect gas over ^ 



kpc scales from gravitational instabilities and spiral arms are longer than local dynamical times, our 
models may nevertheless provide a good first approximation to the effects of star formation feedback 
on local regions within larger gas accumulations. In addition, initial tests we have conducted which 
include magnetic fields (permitting Parker instability) show similar behavior to our unmagnetized 
models. 



As well as producing ~ kpc-scale overdensities, both gravitational instabilities (e.g. Wada et 


aLJ|2002 


Kim et al.||2003 |Kim & Ostriker||2007 ; Agertz et al.||2009 


Aumer et al.||2010t |Bournaud et 


al. 2010 


and spiral shocks (e.g. Kim &; Ostriker 


2006; Kim et al. 


2006, 2010 Dobbs et al. 


2006), 


together with magnetorotational instabilities (e.g. 


Kim et al.||2003 


Piontek k Ostriker||2004 


2005 



2007), drive turbulence in the ISM. In particular, turbulence levels ^ 10 km s~^ can be produced 
by large-scale gravitational instability, and may be important during the highly-transient early 
evolution of disk galaxies. Several of the above numerical models have shown, however, that unless 
energy (representing feedback) is locally injected into massive, high-density clumps that form, the 
result is irreversible gravitational collapse and star formation far exceeding observed rates. Stellar 
feedback therefore appears to be crucial for disrupting bound clouds (thus limiting star formation) 
and maintaining - over many galactic orbits - turbulent ISM levels similar to those observed in 
nearby galaxies. 

The plan of this paper is as follows. In Section 2, we begin by summarizing the theory 
developed in OMLIO and OSll. Section 3 describes the numerical methods and parameters used 
for our time-dependent simulations, and Section 4 presents our model results. These include time 
averages of star formation rates, thermal and turbulent pressures, gas layer scale-heights, thermal 
and turbulent velocity dispersions, and mass fractions of gas components. In Section 5, we use our 
numerical results to test the validity of the physical assumptions and adopted parameters in the 
OMLIO theory. Here, we also demonstrate the balance between turbulent driving and dissipation (as 
in OSll), and quantify the feedback yield relations between P^h and -Pturbj and Ssfr- We compare 
our numerical results to several simple prescriptions for star formation in Section 6. Section 7 
summarizes and discusses our main results. 
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2. Summary of Thermal/Dynamical Equilibrium Model 

In this section, we briefly summarize tlie OMLIO tiiermal/dynamical equilibrium model, high- 
lighting the fundamental assumptions and predictions that we shall test in this work. We then 
draw on OS 11 to outline additional predictions related to the dynamical state and star formation 
rate in disks dominated by turbulent, diffuse gas, and describe how these hypotheses will be tested. 

OMLIO considered a multiphase, turbulent galactic ISM disk with thermal properties mediated 
by stellar heating. The gaseous disk, with total surface density of neutral gas S, is immersed within 
the stellar disk and dark matter halo, whose combined midplane density is given by ps + Pdm = PsA- 
The neutral gas disk is composed of two components: diffuse gas, with surface density averaged 
over large scales Sdiff! and gravitationally bound clouds (GBCs) with surface density averaged 
over large scales (i.e. many individual bound clouds) Sgbc = S — S^iff. The diffuse component 
includes both warm, rarefied gas and cold, dense gas in cloudlets that are not massive enough to 
be gravitationally bound. Star formation takes place within the gravitationally-bound component. 

The first key assumption of OMLIO is that the volume-filling diffuse ISM disk is in force balance 
in the vertical direction. The combined inward gravitational force of the stars, dark matter, and 
gas (both diffuse and GBC components) must be matched by the outward pressure forces within 
the diffuse gas. Averaging the vertical component of the momentum equation over time and in the 
horizontal direction, OMLIO showed that in a state of dynamical equilibrium, Ptot = ^tot de for 



tot.DE 



1 + 2^ + 



iff 



1 + 2 



Sgbc \ ^ 32Qcl,f ^a psd 
SdifT J ttG 



^diff 



1/2- 



(1) 



that is, the total effective midplane pressur^ Ptot must support the weight of the overlying diffuse 
gas in the total gravitational field. Although we use the symbol -Ptot,DE to denote the vertical 
weight, it is important to note that the weight and effective pressure balance only if equilibrium 
holds, and only in an averaged sense. 

In equation ([T]), a is the ratio of total (effective) pressure to thermal pressure in the diffuse 
medium, Q is a dimensionless parameter characterizing the gas density profile {(d = I/tt for a 



,diff/c«; for 



Gaussian profile), = {kT^/ pY^"^ is the thermal speed of the warm gas, and fw 
vth^diR the mass- weighted thermal velocity dispersion in the diffuse gas. The quantity is also 
equal to Pw/ Po for p^ the warm medium density and po the volume-averaged density of the diffuse 
medium (including cold cloudlets, assumed to be in pressure equilibrium with the warm medium) 
at the disk midplane. The mass fraction of the warm medium in the diffuse gas is comparable to 



^As discussed in OMLIO, Ptot is actually a pressure difference between the midplane and the top of the neutral 
layer. Thus, if the cosmic-ray and magnetic scale heights far exceed that of the neutral gas, there is not a significant 
contribution to Ptot from magnetic or cosmic-ray terms (even if their midplane pressures are large), and the weight 
of the diffuse neutral layer must be supported primarily by turbulent and thermal pressure. 



fu, (see OMLIO). In a state of dynamical equilibrium, the midplane diffuse-gas thermal pressure 
^th = Po^^th.diff is equal to 

_ -ftot.DE 



th,DE 



a 



(2) 



(see equation 11 of OMLIO). If the dominant contributions to the total effective pressure are 



(^th,diff + ^z,diff 



)/'^th,difr 



4iff 



/^th.diff 



thermal and turbulent terms with Pturb = Pov'^diS^ then a 
for v^^fiis the turbulent vertical velocity dispersion and cr^^dis the total vertical velocity dispersion 
in the diffuse gas {(Tz,diS is a direct observable for a face-on disk). Note that in equation ([l]), the 
product c^fwCi = Ptot/po, which is equal to cr'^diS turbulent and thermal terms dominate the 
effective pressure, i.e. Ptot ~ Pth + -fturb- 

Next, OMLIO assumed that the diffuse ISM is in a state of thermal equilibrium, in which cold 
and warm atomic phases coexist at a midplane thermal pressure -Pth,TE- In order for the diffuse 
gas to be in the two-phase regime, -Pth,TE rnust fall between the minimum pressure Pmin for the 
cold phase and the maximum pressure Pmax for the warm phase (cf. [Field et al.|[T969| ). Both P^ax 
and Pmin depend linearly on the local radiative heating rate per particle, L, which itself depends 
approximately linearly on the locally-averaged star formation rate surface density, SsfR; if young 
massive stars are responsible for most of the heating. Motivated by detailed modeling of heating 
and cooling in the Solar neighborhood ( Wolfire et al.|[2003 ) and numerical simulations of turbulent 



multiphase gas (Piontek & Ostriker 2005 2007), OMLIO assumed that Pth,TE is comparable to 
the geometric-mean pressure Ptwo = (PminPmax)^/^- Based on the results of |Wolfire et al.| ( |2003[ ), 
OMLIO adopted a geometric mean "two- phase" pressure given by 



Ptwo/^B = 3 X 10=^ cm-3 K 



4G' 



l + 3Z^(S/So)0-4' 



(3) 



Jf\jv/Jf\jv,o is the mean FUV intensity relative to 

2 



where k-Q is the Boltzmann constant, Gq 

the Solar neighborhood value Jfuv,o = 2.2 x 10"'* erg s^^ cm^'^sr^^, Sq = 10 Mq pc~'^ is the 
surface density of neutral gas at the Solar circle (Dickey & Lockman 1990 Kalberla &: Kerp||2009 ), 
and Z'^ is the dust abundance relative to Solar neighborhood value. In the Solar neighborhood, 
Pt^o/kB = 3000 cm-3 K for the OMLIO prescription. 

In a state of simultaneous thermal and dynamical equilibrium, heating and cooling are in 
balance so that Pth = Pth,TE ~ Ptwoi and vertical forces are in balance so that Pth = Pth.DE- 
With Ptwo oc G'q oc Jpuv 5]sFR, the surface density of star formation should be proportional 
to Pth- Thus, equating ^ and ([s]) yields an expression for the star formation rate, with Esfr 
proportional to the right-hand side of equation ([T]) - i.e. to the weight of the diffuse gas layer in the 
total gravitational field. In low-density outer-disk regions where the diffuse gas dominates GBCs 
(Sdiff — >• S and Sqbc — 0), an approximate form for Ssfr is then given by 
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1/2 , , 1/2" 



aVlOMopc-^y \ a / VO-1 pc 



(4) 



(see eqs. 22 and A13 in OMLIO). The numerical coefficient in e quation Q i s calibrated based on 
the local Milky Way value Ssfr.o = 2.5 x lO^^ Mq kpc"^ yr"! ( |Fuchs et al.|2009[ ). 



In the case when Sgbc/^^ is non- negligible, in order to obtain a closed set of equations, OMLIO 
made the additional assumption that star formation within GBCs has a gas consumption timescale 
isF,GBC so that 

SSFR = ^^=^^^. (5) 
iSF,GBC *SF,GBC 

If GBCs have relatively uniform properties, then isF.CBC will be relatively constant. By equating 
([2]) and ([3]), and combining with equation ([5]), OMLIO obtained a cubic equation that can be solved 
for SgFR as a function of S and psd in the general case; an approximate form is given by 

-1 



^SFR 



(6) 



S SsFR,low 

(see eqs. 23 and A14 in OMLIO). Note that for low surface density outer disks, equation Q is 
recovered and Ssfr is independent of tsF.GBC ~ i-e. the star formation rate becomes independent 
of the rate at which gas in GBCs collapses to make stars. 



OMLIO took tsF.GBC = 2 Gyr based on the empirical linear correlation (Bigiel et al. 2008) 
between the molecular mass in CO and the SFR at 750 pc scale for a set of disk galaxies (at moderate 
E ^ 100 M0 pc~^), and adopted a w 5 and fw ~ 0.5 as typical values based on observations of 
the Milky Way and other well-studied disk galaxies. If the same set of parameters is adopted for 
all galaxies (note that the dependence on fw/a in equation [i] is weak: Ssfr oc (/^/q;)°'^), Ssfr is 
a function of just S and ps^- OMLIO applied this formulation to azimuthally-averaged data for a 
sample of spiral galaxies to predict SgpR as a function of galactocentric radius R. The resulting 
predicted profiles of Ssfr are overall in remarkably good agreement with the observations. For a 
few galaxies, however, observed values of SgpR are offset from the prediction. The difference may 
owe to different values of a, fw, and/or tsF,GBC from the adopted values, or to effects associated 
with azimuthal averaging when there is strong spiral structure. It should also be noted that there 
are still significant uncertainties in the observations, which might lead to offsets with respect to the 
theory. Empirical determinations of S and tsF.GBC are uncertain since some gas may be undetected 
in both 21 cm and CO lines, and since the conversion factor Xqo from CO to H2 can vary by a 
factor ~ 2 {Xqo varies even more at low metallicity, and where S ^ 100 Mq pc~^). The age of 
the young-star population, as well as the treatment of extended vs. concentrated tracers of star 
formation, can also affect the empirical estimates of SgpR. In addition, as discussed by OMLIO, 
values of psd are uncertain as stellar disk thickness estimates for face-on galaxies are obtained via 
scaling relations rather than being directly measured. 

In this paper, we focus on the low-S case, corresponding to outer disks where the gas is 
primarily diffuse and atomic. In this regime, SsFR is predicted to depend on a and fw but not 
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on tsF,GBCj according to equation ([4j). Using our numerical simulations, in which S and psd are 
independent variables, we can directly test the primary assumptions of the OMLIO theory. Since 
we can measure a, fw, ^dis (and Sqbc = S — Sdifi) together with Pth from the simulation outputs 
for any model, we can test whether the measured midplane thermal pressure in fact agrees with 
the dynamical equilibrium value Pth,DE predicted by equation Q . We can also investigate whether 
the measured midplane Pth is close to Ptwo following the hypothesis of OMLIO that the system 
evolves to a state of thermal equilibrium having both a warm and cold atomic phase. Similarly, we 
can test whether the sum of the measured thermal and turbulent pressures Pth + Po^'^diS ~ ^tot is 
consistent with the dynamical equilibrium prediction of equation ([T| (since the present simulations 
do not include magnetic fields, cosmic rays, or radiation pressure, these terms do not enter Ptot)- 
Further, we can check whether our numerical results for a and agree with empirically-estimated 
values, and explore how much variation in a and fw there is among models with different S and 
Psd- Finally, we can compare the value of Ssfr from the simulations with the theoretical prediction 
based on simultaneous thermal and dynamical equilibrium (cf. equation |4]) . 

In addition to testing the theory of OMLIO, we can use our numerical simulations to test more 
general ideas related to the self-regulation of star formation, as introduced by OSll. We consider 
the situation in which the ISM is dominated by diffuse gas, so that ScBc/^^difr — ^ and Sdis — )• S. 
We also assume the effective pressure is dominated by thermal and turbulent term^ and take 



Cd ~ I/tt and cif^a = cj^ 



,difr 



aj so that equation ( 1 ) for the weight becomes 



Pi 



t:GT? 



tot,DE 



1 + 



1 + 



1/2- 



(7) 



A simplified expression for Ptot,DE) within 20% of equation ([7]), is 



tot,DE 



+ E(2Gpsd)'/V, 



lO^/cB cm^'^K 



0.33 



10 M0 pc- 



+ L4 



Psd 



(8) 



1/2 



O.lMopc^^y VlOkms 



0-2 



10 Mo pc-2. 
The vertical dynamical equilibrium equation is 

-fth + -fturb = PO'Wth.difr + PO^^idiff = P^^'^l = -ftot.DE- 



,-1 



(9) 



As noted above, it is expected that Pth = Pa'^th diff ^sfr in a state of thermal equilibrium. 
In addition, OSll argued that if mechanical feedback from star formation provides the dominant 



^That is, we assume cosmic ray, magnetic field, and radiation effects are unimportant - see OMLIO and OSll for 
an evaluation and discussion of these. 
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contribution to the vertical turbulent motions, then the turbulent pressure Pturb should also scale 
roughly linearly with Ssfr, as 

-fturb = /p^^5]sFR- (10) 

Here, is the mean radial momentum injected by each massive star, m* is the total mass in stars 
formed per massive star, and the order-unity coefficient fp parameterizes the details of turbulent 
momentum injection and dissipation. When turbulence dominates the pressure and self-gravity 
dominates the vertical weight, equations ([7]), ([9]) and (10) with /p 1 combine to yield a prediction 



that SsFR ~ 27rGS^m*/p*. OSll found that this prediction is in good agreement with both 
numerical simulations (for a cold-gas dominated ISM) and with observations of molecule-dominated 
starburst regions with S ^ 100 Mq pc~^. 

More generally, if star formation is responsible for both heating and driving vertical motions in 
the diffuse ISM, we expect the thermal and turbulent pressure contributions to scale roughly linearly 
with SgpR. Normalizing relative to convenient dimensional units for observational comparison, we 
can define 

~ ^th—z^T^^ =2 ZT (11) 



103 cm-3 K 10"3 M0 kpc-2 yr- 

-Pturb/^B ^ ^SFR 

103cm-3K - ''*"'-'^lO-3M0kpc-2yr-i' ^ ^ 

The parameters ryth and ?7turb are yield coefficients that measure the efficacy of feedback. For 
the fiducial parameters adopted in OMLIO, r?th = 1.2[0.25 + 0.75Z^(S/10 Mq pc-2)0.4]-i^ ^j^ere 
the factor in square brackets is unity in the Solar neighborhood. For the fiducial value = 
3000 km s-i adopted in OSll (assuming supernovae are the most important sources of momentum), 
??turb = 3.6/p. Note that with the heating and turbulent driving yield coefficients as defined in 



equations (11) and (12), a = (Pth + -Pturb)/-fth = 1 + ??turb/^th if only thermal and turbulent 
stresses contribute to the effective midplane pressure. We thus expect J^th+^turb ~ 1.2+3.6 ~ 5 and 
a ~ 1 + (3.6/1.2) = 4 under conditions similar to the Solar neighborhood. The latter is comparable 
to the value a = 5 adopted in OMLIO for comparisons of equation ^ with observations of Ssfr- 
By exploring the relations between the measured values of Pth, -Fturb) and Ssfr in our simulations, 
we can numerically evaluate r^th and ?/turb) testing whether these quantities (and therefore q) are 
indeed near-constant. 

Combining equations ([9]), (11) and (12), the self-regulated star formation rate in a diffuse-gas- 
dominated region where the pressure is controlled by energy and momentum feedback from massive 
stars has the form 

EsFR = 2 X 10-3 Mq kpc- yr- (^^Ml^^ " W^^' (^3) 

For outer-disk regions, equation ([7]) or ([s]) may be used for the ISM weight Ptot,DE- In galactic- 
center regions where the bulge potential exceeds that of the disk, p^d Ph/'i for p}, the bulge stellar 
density (see OSll). 
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For very dust-poor systems, FUV radiation escapes more easily from star-forming regions and 
penetrates further in the diffuse ISM, which may make the heating yield r/th comparable to or even 



larger than r/turb (see OMLIO and Bolatto et al. 2011 ). Alternatively, in regions where S is extremely 
high and reprocessed IR radiation is trapped, radiation pressure becomes important and a term 



r/rad ^ Skir would be included in equation (13). Since the cosmic ray and magnetic pressures 



presumably increase with higher Ssfr in analogy with equations (11) and (12), corresponding 



feedback terms could be included in equation (13), with the values of ?7cr and r/mag appropriately 
taking account of differing vertical scale heights compared to the neutral, star-forming gas (see 
OSll). 

Using our present simulations, we can test whether the generalized feedback- regulated star 
formation prediction SsFR oc -Ptot,DE is satisfied. We will also compare our results to the power-law 
form SgFR oc traditionally used in fitting observations, and to the form SgpR ^Po"^ that 

is frequently adopted in numerical simulations of galaxy formation/evolution in the cosmological 
context. 



3. Numerical Methods and Models 

3.1. Basic Equations 

The numerical models of this paper investigate thermal and dynamical evolution of gas in 
a vertically stratified, differentially rotating, self-gravitating galactic disk under the influence of 
interstellar cooling, heating, and radiative and mechanical feedback from star formation. We set 
up a local Cartesian frame whose center is located at a galactocentric radius Rq and rotates with 
an angular velocity = ^}{Ro). In this local frame, x = R — Rq, y = Roi4' ~ aiid z represent 
the radial, azimuthal, and vertical coordinates, respectively. Our simulation domain is a two- 
dimensional rectangular region with size x in the x - z plane with ?/ = (hereafter XZ 
plane), representing a radial- vertical slice of the disk, although we implicitly consider the thickness 
-Ly(<C Lx,Lz) in the y-direction for the purposes of computing star formation rates and momentum 



feedback (see Section 3.2.1 ). We include nonzero velocity in the y-direction in order to treat epicyclic 
motions self-consistently. The equilibrium background velocity relative to the center {x = z = 0) 
of the simulation domain is given by vq = —qQxy, where q = —{dln^l/dlnR)\f{^ is the local 
dimensionless shear rate. In terms of q, the epicycle frequency k is given by = (4 — 2q)Q'^. We 
assume a flat rotation curve so that q = 1 and k = ^/2^}. 

We expand the basic equations of hydrodynamics in the local frame, neglecting terms arising 
from the curvilinear geometry. The resulting shearing-sheet equations (e.g. Kim et al.|2002 Piontek 



fc Ostriker|[2()07| ) are 

1^ + V • (pv) = 0, (14) 
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dv 



+ V • Vv 



1. 



-VP - 211 X V + 2gfi\x - V$ + gsd, 



(15) 



de 



at 



+ v 



(ev) = -PV • V - p£ + /CV^T, (16) 

= 47rGp, (17) 

where is the self-gravitational potential of the gas, gsd is the external gravity from the stellar 
disk and the dark matter halo, pC is the net cooling function, and fC is the thermal conductivity. 
Assuming that the gas is predominantly atomic and has cosmic abundances, P = l.lnk-QT is the 
gas pressure where n = p/{lAmp) is the number density of hydrogen nuclei. We adopt an ideal 
gas law so that the internal energy density is given by e = P/(7 — 1) with index 7 = 5/3. For the 
external gravity, we take the simple form 



gsd = -47rG/9sd^;z, 



(18) 



where /?sd is the midplane density of the stellar disk plus that of the dark matter halo. Since the 
scale height of the gas is much smaller than those of the stellar disk and the dark matter halo, gsd 



given in equation (18), corresponding to vertically- uniform p^d, is a reasonable approximation in 



studying dynamics of the gas. 

The net cooling function per volume is given by pC = n[nA(T) — F]. For the cooling rate of 



the diffuse ISM, we adopt the fitting formula obtained by Koyama &: Inutsuka (2002): 

-1.184 X 10^ 



A(r) = 2 X 10~^^exp 



T + 1000 



+ 2.8 X 10 



-28 



Texp 



-92 



erg cm s 



(19) 



with temperature T in degrees Kelvin. Cooling at low T is dominated by the 158/xm fine-structure 
line of C II, whereas cooling at high T is dominated by Lya line emission; both lines are collisionally 
excited. The heating rate F is dominated by the photoelectric effect on small dust grains and 
polycyclic aromatic hydrocarbons (PAHs) by FUV photons with energy 6 eV < hv < 13.6 eV 
( Bakes &: Tielens|1994 ). The diffuse FUV radiation field, with intensity JfuV) is produced by young 
O and B stars and therefore should depend on recent star formation. We thus allow F to vary with 
time, while keeping F uniform throughout the simulation box (i.e. Jpuv is treated as spatially 



constant). We follow Koyama k, Inutsuka (2002) in adopting a fiducial heating rate in the Solar 
neighborhood Fq = 2 x 10~^^ erg s~^. In thermal equilibrium [pC = 0) for this cooling function, two 
stable phases co-exist for a range of densities and pressures: the maximum pressure for the warm 
phase is Pmax/A^B = 5.5 x IQ^{T/Tq) cmT^ K occurring at T^ax = 5000 K and ni = 1.0(F/Fo) cm'^, 
and the minimum pressure for the cold phase is Pmin/^B = 1-8 x 10'^(F/Fo) cm~'^ K at T^\a = 188 K 
and 712 = 8.7(F/Fo) cm~^. The two-phase pressure is thus given by Ptwo/^B = (PminPmax)^''^/fcB = 
3.1 X 10'^(F/Fo) cm~^ K. For Solar- neighborhood conditions, Ptwo is essentially the same as adopted 
in OMLIO, Ptwo/^B = 3000 cm^'^ K (see equation We describe our prescription for connecting 



F with the (time-dependent) star formation rate, including metagalactic FUV radiation, in § 3.2.2| 



Thermal conduction plays an important role in the development of thermal instability (TI). 



Conduction not only sets the critical wavelength (the "Field length") of TI (Field 1965), but also 
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determines the thickness of interface layers between cold and warm phases (Begelman & McKee 



1990 ) . Inclusion of thermal conductivity is therefore essential to resolve TI in numerical simulations 
dKoyama k Inutsuka||2004t |Piontek Ostriker||2004t |Kim et al.||2008[ ). A realistic value of thermal 
conductivity in the diffuse ISM at T < 10^ K is /C ~ 2.5 x IO^tVs erg s'^ cm"^ K"^ ( |Parker 
1953). The corresponding Field length is then Xp ~ 0.2 pc for the typical density n = 1 cm~^ 
and temperature T = 10^ K of the thermally unstable gas, which would require an extremely fine 
numerical grid Ax ^ Air/3 in order for TI to be resolved. In addition, hydrodynamic simulations 
involving supersonic turbulence inherently suffer from a significant level of numerical diffusion 



Gazol et al. 


2005 


Kim et al. 


2008 



Ax is extremely small. Adopting a realistic value of /C is therefore prohibitively expensive for 
multi-dimensional simulations in kpc-scale numerical boxes. Fortunately, however, dynamics on 
larger scales are not sensitive to the exact conduction scale, similar to large-scale dynamics in 
supersonic fiows being insensitive to the exact thickness of shocks. In this paper, we therefore 
adopt a numerical conductivity of /C = 4 x 10'' erg s~^ cm~^ K~^/[l -|- (0.05 cm~^/n)] as in Koyama 



&: Ostriker (2009a), which enables us to resolve the Field length numerically, and limits thermal 



conduction in low-density regions. 



We solve the time-dependent partial differential equations (14)-(17) using a modified version 
of the Athena code (Stone et al. 2008: Stone &; Gardiner 2009| ). Athena employs a single-step, 
directionally unsplit Godunov method for (magneto) hydrodynamics in multispatial dimensions, 
providing several schemes for integration in time, spatial reconstruction, and solution of the Rie- 



mann problem. We use the van Leer algorithm (Stone &: Gardiner 2009) for integration, with 



piecewise linear reconstruction and the HLLC Riemann solver. We solve the net cooling function 



based on implicit time integration using Simpson's rule (e.g., Koyama & Ostriker 2009a) with a 
limit for the maximum temperature change of 50%. We also use an explicit conduction solver for 
isotropic thermal conduction, and revert to first order flux updates if a negative density appears 
during the higher-order update ( Lemaster &: Stoiie||2009 ) . The gravitational potential is calculated 
using fast Fourier transforms in disk geometry with vacuum boundary conditions in the z-direction 
( Koyama &: Ostriker|2009a ) . At the x-boundaries, we apply shearing-periodic boundary conditions 
(Hawley et al. 1995). In the z-direction, we adopt periodic boundary conditions for the hydro- 



dynamic variables so as to maintain a constant mass within the domain. By running comparison 
models using outflow boundary conditions in z, we have checked that the boundary conditions do 
not affect the simulation outcomes significantly. 



3.2. Prescription for Star Formation Feedback 

In our simulations, self-gravitational collapse and ensuing feedback from star formation control 
both thermal and dynamical evolution of the model ISM. We consider both mechanical (momentum 
input) and radiative (thermal energy input) feedback effects. Mechanical feedback drives turbulence 
that supports the disk in the vertical direction, while radiative feedback affects the thermal pressure 
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by changing the heating rate. In this section, we detail our prescription for star formation feedback]^ 



3.2.1. Mechanical Feedback 



Star formation in our models occurs only inside clouds where the gas density is larger than 
a critical value. The threshold density pcr should be large enough for star formation to occur 
only in self-gravitating regions. In addition, these self-gravitating regions should be resolved on 
the grid, i.e. the Jeans wavelength Aj(/?cr) = [''^c^r/ {Gpcr)]^^"^ should exceed the grid spacing Aj; 
(taken to be 1 pc in our models), where Ccr denotes the thermal speed at the threshold temperature 
TcT- Since the cooling time is very short, dense clouds are generally in thermal equilibrium, and 
Tier = r/A(Tcr). Equation (19) then yields 



Aj ^ i.4r3/4e-46/^-(r/ro)-i/2 pc. 



(20) 



for T < 100 K. For a fixed Aj, we obtain T^r (and hence n, 
power-law fit for A j = 2.7 pc gives ricr ~ 500(r/ro)'^'^ 
density for star formation in our simulations. Although slightly lower threshold density would be 



CT) as a function of T /Vq. A simple 
which we take as the threshold 



needed to meet the Truelove criterion Aj/Ax > 4 (Truelove et al. 1997, 1998) and limit artificial 



fragmentation in collapsing clouds, our choice is acceptable in the current context since our aim 
is not to follow cloud collapse and fragmentation but instead to disperse self-gravitating clouds by 
turning on star formation feedback, as explained below. 

Not all clouds with p > pcr immediately undergo gravitational collapse and star formation, 
since the star formation efficiency and the computational time step should be considered as well. 
Let us consider a star- forming region with density p > pcr- Assuming that our simulation domain 
represents a two-dimensional slab with thickness Ly in the y-direction, the mass in the cloud above 
the threshold is Md = Ly Jp^p pdxdz. For the thickness of the slab, we take Ly = 2rsh, where r^^ 
is the initial radius of an SN shell explained below. This choice of Ly is due to the fact that the 
most significant feedback in the simulation domain comes from SN events occurring within 2rsh in 
the y-direction. The SFR expected from the cloud is 



eg 



%(P) 



(21) 



where eg is the star formation efficiency per free-faU time, %(p) = {Zt: / {?,2Gp)Y'^ . We take 
eg = 0.01 as a fiducial value consistent with theory and observations (Krumholz & McKee 2005 



''other recent numerical studies of the ISM have used somewhat different prescriptions for radiative and mechanical 
feedback from those we adopt. For example, Joung et al. (20091 adopted V oc Eg^^ together with type-II SN rates 
scaling as Esn oc Y,], 



Agertz et al. ( 2009 1 included feedback from supernovae based on a volumetric star formation 



-■gas I 

rate psfr oc Pg^s but did not include diffuse UV heating; and 
declines exponentially outward, but did not include mechanical feedback from supernovae. 



Tasker 



(2011 1 adopted a photoelectric heating rate that 
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Krumholz Sz Tan||2007 ). The probable number of massive stars to form within the cloud in a time 
interval At is then given by 

""^At, (22) 



where is the total mass of stars in all masses formed per massive star. We define massive stars 
as those that undergo supernovae, and adopt m=K = 100 for all simulations consistent with the 



initial mass function of Kroupa (2001). For a given computational time step At, J\f* calculated 
from equation (22) is typically ~ 10~^ — 10"'^ (as small as ~ 10~^ immediately after SN explosions 
due to small time step), much smaller than unity. Therefore, in zones where p > pcr we generate a 
uniform random number M G [0, 1) at each time step, and turn on feedback only provided M^^ > M. 

We implement mechanical feedback from star formation in a very simple way, by injecting 
momentum in the form of an expanding spherical velocity distribution to represent the radiative 



stage of a SN (cf., Shetty &: Ostriker 2008). As the initial radius of the shell in three dimensions, 
we take = 10 pc, corresponding to the SN shock radius at the shell formation time (Cioffi et 



al. 


1988 


Koo Sz Kang 


2004 



randomly in the range \yos\ < ''sh) so that the initial shell radius in the XZ plane (at y = 0) is 



'sh 



1 /2 

y^g) , varying between and rgh- We use a random number to choose the value of 



UoS for each feedback event. When a feedback event occurs, we first redistribute mass, momentum, 
and thermal energy within a circular region of radius -Rmax by taking spatial averages. We then 
add to the momentum density in the x- and z-directions according to 



PVsh,2D 



Pmax 

0, 



R > Rrnax) 



(23) 



where R is the position vector with respect to the center of the SN sphere in the XZ plane, and Pmax 
is the momentum density at R = i?max- By requiring the mean momentum input from equation 



(23) (averaged over i/os) is equal to the outward momentum that a three-dimensional shell would 



have, one obtains Pmax = 15p*/(32rgjj), where is the total radial momentum in three dimensions. 
In all simulations, we take = 3 x 10^ Mq km corresponding to the late stages of a single 
SN with energy Esn = 10^^ erg (Cioffi et al. 1988). The velocity profile v{R) oc R^ is chosen to 
guarantee an initially divergence- free velocity field at R = 0. 

We note a few caveats that should be kept in mind regarding our simplified prescription for 
star formation feedback. First, as our main focus is on the diffuse gas component (which dominates 
by mass), our treatment does not attempt to follow the evolution and destruction of star- forming 
clouds in detail. Thus, we do not introduce a time delay prior to the momentum injection, or 
separately model effects of expanding H II regions or winds (the former was previously considered 
in |Koyama &: Qstriker|2009a[ which found that only relatively low levels of turbulence were induced 
in the diffuse ISM). In this first study, our goal is primarily to incorporate turbulent driving in 
the diffuse ISM at a realistic level for a given star formation rate, which is accomplished by simply 
injecting momentum impulsively. Future work should improve this treatment, but experience with 
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numerical models of turbulent giant molecular clouds has shown that much astronomical insight 



can be gained even when idealized treatments of turbulent driving are adopted ( Mac Low & Klessen 



2004 McKee & Ostriker 2007). 



Second, although our feedback treatment aims to model turbulent driving in the neutral 
warm/cold ISM that is induced by SNe, our approach does not attempt to model the high- 
temperature interiors of SN remnants themselves. Previous work has shown that it is difficult 
to model SN explosions by injecting thermal energy in large scale simulations because of overcool- 
ing: radiative energy losses are too rapid due to lack of spatial resolution ( Katz|1992 ). At resolution 
levels that are affordable, far too little thermal energy ends up being converted to kinetic energy; 
instead it is radiated away. In order to avoid overcooling, in some simulations radiative cooling is 



artificially turned off until blast waves have developed (e.g., Thacker Sz Couchman 2001; Agertz 



et al. 2011), or the initial sizes of regions where SN energy is injected are set such that the gas 



temperature T ~ 10^ K, where a dip is present in the cooling function (Joung &: Mac Low 2006). 



For simulations such as ours which include self-gravity, SN events occur within very dense regions. 
Since the cooling rate is proportional to the square of the gas density, experiments we conducted 
with thermal energy injection and a coronal-gas cooling function showed that the cooling time was 
still unrealistically short at the resolution of our simulations, even if we adjusted the gas tempera- 
ture to the dip of cooling function. Thus, although hot gas created in SNe may be quite important 
in many ways (including driving galactic winds), the present models focus just on the warm/cold 
ISM and star formation, and leave the interesting issues of the hot ISM for future work. 



3.2.2. Radiative Feedback 

Since the photoelectric heating rate is proportional to the intensity of the FUV radiation field, 
we simply take T oc JfuVi with a proportionality constant depending on the heating efficiency 



of small grains and PAHs (see e.g., Bakes & Tielens 1994). There are two sources of the FUV 



radiation field in outer disk: JFUV.iocah the FUV radiation emitted by recently- formed OB stars 
locally in the disk, and Jpuv.meta, the metagalactic FUV radiation field. Radiation originating in 
the inner regions of the galaxy could also reach the outer galaxy, but this contribution is smaller 
than the local radiation unless the optical depth is very low. 

If FUV escapes into the diffuse ISM from star-forming regions at the midplane at rate per 
unit area SfuVj then Jpuv = 5]fuv[1 — E2{t±/2)]/{4:-7tt±) for t± = Skfuv the optical depth 
through the diffuse neutral ISM, and E2 the second exponential integral. As the radiative transfer 
factor depends only logarithmically on 1/t± at low optical depth, for simplicity OMLIO adopted 
Jfuv 5]fuv oc Ssfr for application to galaxies with dust abundance not far from Solar and a 
moderate range of diffuse-H I surface densities. In galaxies with very low dust abundance, UV may 
escape much more easily from star forming regions, and also travel further through the diffuse ISM. 
This would lead to an increase in both Sfuv/^sfr and Jfuv/^fuv relative to the Milky Way, so 



that the ratio Jfuv/^^sfr could be much higher than in the Solar neighborhood. Bolatto et al 
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(2011) found that the warm H I and star formation content of the SMC indeed appears to require 



a higher ratio of Jfuv/^^sfr than in normal disks hke the Milky Way. 

In this work, we assume the heating rate due to local FUV scales with the local star formation 
rate as T/Tq = /rad^^SFR/^^SFR.O; where Ssfr,o = 2.5 x 10~'^ Mq kpc~^ yr~^ is the SFR surface 
density in the Solar neighborhood (Fuchs et al. 2009) and Fq = 2 x 10~^^ erg s~^ (Koyama & 



Inutsuka 2002). The parameter /rad thus implicitly includes the normalized heating efficiency 



of the FUV radiation, allows for additional forms of heating such as X-rays (see Wolfire et al. 
1995, 2003), and would vary depending on details of radiative transfer. Note that /rad = 4/[l + 



3(Z^E/10 Mq pc-2)0-4] is adopted in OMLIO based on the fit in [Wolfire et al.| ( |2003[ ); this has 
/rad = 1 in the Solar neighborhood. 



The total volumetric heating rate is then written as 



/ 



rad 



SsFR \ _|_ / '^^FUV^ 
5^SFR,0/ \ -'fUV,! 



mcta 



(24) 



Note that the heating by the metagalactic FUV given by the second term in equation ( 24 ) provides 
a minimum heating rate when Ssfr is extremely small. We adopt Jfuv meta = 0.0024 Jfuv o 



(Sternberg et al. 2002), so that in practice Jpuv.meta is negligible in most cases. The cooling and 



heating rates we adopt give geometric-mean two-phase pressure equal to 

Ssfr 



-Ptwo/ h 



1.2 X 10^ cm^^K/ 



rad 



(25) 



.10-3 Mq kpc" 

Thus, comparing to equation (11), if we were to find Pth = -ftwo for the mean midplane thermal 
pressure, it would imply r/th = 1-2/rad for the dimensionless heating-feedback yield coefficient. As 
we shall show in Section 5.2, Pth at the midplane is in fact between +10% and —40% of Ptwo, so 



that remains very close to 1 x /^ad- 

In order to change the heating rate self-consistently, we need to calculate the recent SFR at 
each time step. We do this by counting the number of the recent SN events, so that the SFR surface 
density is calculated by 



SSFR 



(26) 



where tbin is the time bin over which the SFR is averaged, and Ng-^ denotes the total number 
of SN events that occurred during the time span {t — tbin, t). We note that Ssfr corresponds 
to a space and time average of divided by the surface area. Since only recent star formation 
contributes to gas heating via FUV radiation, if the simulations were in three dimensions and 
optical depth effects were included, the averages should be taken at least over tFUV x i'^d?) to cover 
the whole domain of influence, where fpuv ~ 10 Myr is the FUV luminosity-weighted lifetime 
of OB stars (Parravano et al.|2003) and d ~ 200 pc/(S/10 Mq pc~^) is the effective in-plane 



distance for radiation to travel^ However, our simulation domain represents a radial- vertical slab 



^ The effective in-plane distance for FUV radiation to travel is given by d ~ 2///(Skpuv) where H is the scale 
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with effective thickness Ly = 2rsh = 20 pc in the y-direction, with Ly <ti d. Since the size of 
our domain in the x-direction is large enough (L^ ^ d), it is desirable to take a temporal bin at 
least ibin ~ t¥w{d/ Ly) ~ 10tFUv/(S/10 M© pc~^) in order to limit stochasticity in the heating 
rate. We thus set ibin equal to a half of the orbital period (see below for definition). Since our 
set of model parameters is chosen to maintain oc S, this implies ^bin oc With this choice, 

tbin(S/10 Mq pc~2) ^ 100 Myr ~ lOtpuv- 



3.3. Model Parameters 

Since the feedback parameters are all specified, we now turn to the disk parameters. Our 
initial conditions for the gaseous disk consist of warm-phase gas with uniform thermal speed — 
7 km s~^. The gravitational susceptibility of the disk depends on three parameters: gas surface 
density S, the angular velocity of galactic rotation 0, and the stellar plus dark matter density at 
the midplane p^d- Both E and enter the Toomre stability parameter 

Qinit = (27) 

while /9sd determines the degree of vertical disk compression induced by the stellar disk and dark 
matter halo. It is convenient to define 

2^;;^ "°-^%10 Mq pc-y V7km s-J VO-05 Mq pc-3 J ' 

which measures the relative strengths (in the vertical direction) of gas self-gravity and the external 
gravity from stars and dark matter ( Kim et al.|2002 >. For Solar-neighborhood conditions, sq ~ 0.3. 



Assuming sq <C 1, the equilibrium density distribution is a Gaussian profile 

p{z) = poeM-zV^H^), (29) 



where po = S/[(27r)i/2i?^] and 

(47rGpsd)^/2 - ^"^^^ V7kms-iy V^-OS M© pc-3 



= . = 134 PC ( — ^ ) ( ] "\ (30) 



is the scale height. 

To simulate disk evolution in a range of environments systematically, we vary S and /?sd 
while keeping Qinit = 2 fixed, so that the angular velocity at the center of the domain varies as 
Q = 28 km s~^ kpc^^(S/10 pc~^). We consider four main series of models: QA, QB, S, and 
G. The model parameters are summarized in Table [l| In Series QA and QB, psd varies as psd oc 



height of the gas disk and kfuv ~ 1 — 2 x 10 cm^( H atom) ^ ~ 0.1 pc^ is the dust opacity in the FUV 

band. By taking H ~ 100 pc, we have d ~ 200 pc/(E/10 Mq pc"^). 
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so that the stellar Toomre parameter Qs oc il/^/p^ oc Q^/sq implicitly has the same value for all 
members of each series. For the QA series, sq = 0.28 and for the QB series sq = 0.07. Thus, models 
in Series QB have four times larger /Osd (i-e. a more confining stellar vertical potential) than those 
with the same S in Series QA. The model Series QA and QB represent conditions typical in disk 
galaxies at different galactocentric radii, from mid-disks (i.e. slightly inside the Solar circle) to far 



outer disks (e.g., Koyama & Ostriker 2009a[)p| For Series S, we fix and vary S to explore the 



effect of the gas surface density independent of the strength of the external vertical gravity. In Series 
G, S and Q are held constant, while psd varies; this allows us to isolate the effect of the external 
vertical gravity. Our fiducial model is Model QAlO with S = 10 Mq pc~^, = 28 km kpc~^, 
and Psd = 0.05 Mq pc~^; this model is similar to the Solar neighborhood. The corresponding 
orbital period is torb = 27r/il = 220 Myr(Q/28 km s"^ kpc"^)~^ = 220 Myr(E/10 M© pc^^)-!^ 
which we use as the time unit in our presentation. 

The model series above all have the same feedback parameters. In addition, we consider Series 
R, in which /rad is varied to explore the effect of varying heating for a given Ssfr. All other 
parameters in Series R are the same as Model QAlO (which has /rad = !)• We ran four models 
labeled R02, R05, R25, and R50 with /rad = 0.25, 0.5, 2.5, and 5.0, respectively. Since the ratio 



of local heating rate to local SFR surface density F/Esfr oc /rad (see equation 24), larger /^ad 
implies a higher heating rate for a given SsfR; corresponding to lower shielding (e.g. from lower 
dust abundance) than in the Solar neighborhood. Smaller /j-ad corresponds to higher shielding. 
In reality, /rad should depend on both dust abundance and the total column of gas, since both of 
these can affect shielding. For the present study, we simply treat /rad as an autonomous variable 
in order to explore effects of varying shielding (or heating efficiency, which for present purposes is 
equivalent). 

For the vertical extent of our simulation boxes, we take = ^Hy^ (this varies depending 
on the model; see Table [T]). In the horizontal direction, we take Lx = 512 pc as the standard 
value. In order to check the effect of the box size, we have run Model QA10x2, which has the same 
parameters as Model QAlO except the horizontal box size is extended to = 1, 024 pc; this model 
confirmed that overall evolution and statistical properties are indeed similar. We vary the number 
of zones from model to model to make the grid spacing Ax = Az = 1 pc for all the models. In 
order to seed TI, isobaric perturbations consisting of a Gaussian random field with flat power for 
1 < kLz/2'K < 8 and zero power for kLz/2TT > 8 are added to the initial density and temperature 
distributions. The amplitude of the initial perturbations is set to 10% of the midplane density. 
We evolve each model until i/torb = 3, well beyond the time required for the system to reach a 
quasi-steady state. 



® Very far outer galaxies with negligible stellar disks and only dark matter contributing to psd oc R ^ would have 
So — 2(7rGE)^(c„,f2dm)^^, which could reach unity, but these conditions are not studied in the current work. 
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3.4. Classification of Gas Components 



Before describing the simulation results, we establish terminology for the various gas compo- 
nents we shall discuss. In the neutral ISM, gas in GBCs and diffuse gas are distinguished based 
on whether the gravitational energy and total pressure significantly exceed that of the surround- 
ing gas at similar z, or not. In general, the GBC component consists of the population of giant 
molecular clouds (GMCs), including both molecular gas inside GMCs and dense atomic shielding 



layers. Observations of the Milky Way (Solomon et al. 1987; Heyer et al. 2009 Roman-Duval et 



al. 2010) and Local Group galaxies (Bolatto et al. 2008) have reported that GMCs have similar 
surface densities Sgmc ~ 100 M© pc~^, corresponding to ncMC ~ 40 cm~^(MGMc/10^ M©)"-^/^. 
Since we do not take into account radiative transfer and formation of hydrogen and CO molecules 
explicitly, we cannot directly identify structures in our models that would be observed as GMCs. 
In this work, we simply define gas with n > uq-qc = 50 cm~^ as being within the GBC component, 
since observed GMCs have comparable densities. We emphasize that this classification is essentially 
a nomenclature shorthand, allowing us to refer to the densest gas as the "GBC component". The 
designation of gas as "GBC" or "diffuse" component is not used in any way within the simulations 



themselves. We note that the density threshold for star formation (see section 3.2.1 ), which is much 
larger than uq^q, ensures that star formation in our numerical models takes place only within the 
GBCs. 

The diffuse component, defined as gas with n < uqbCj consists of thermally-stable cold and 
warm phases as well as a thermally-unstable phase. We classify the phases of the diffuse component 
based on its density rather than temperature such that it is warm gas if n < ni , cold gas if n > n2 , 



and unstable gas if ni < n < ?7-2 (see definitions of rii and n2 following equation 19). Note that ni 
and 712 depend on T (and hence Ssfr) and thus vary with time. In what follows, /gbc and f^is 
denote the mass fractions of GBC and diffuse components in the whole gas, respectively. Similarly, 
the mass fractions of cold, unstable, and warm phases within the diffuse component are represented 
by fc, fu, and respectively. Note that /qbc + /diff = 1 and /c + /« + /«; = 1- 



Simulation Results 



In this section, we describe results of our numerical simulations. Our models evolve in a 



generally similar manner to those of Koyama Sz Ostriker ( 2009a ) , which also included self-gravity 



radiative heating (at fixed F) and cooling, and feedback from star formation. In the models of 



Koyama Sz Ostriker (2009a), only feedback associated with H II regions was considered. H II 



regions were modeled by applying intense heating in dense enough regions that met criteria for star 
formation; expansion of the overpressured gas provided turbulent driving. Since SN explosions are 
more energetic than expanding H II regions, however, our present models achieve a higher (more 



realistic) level of turbulence at saturation than those in Koyama Sz Ostriker (2009a). Also, the 
variable radiative heating rate in the present simulations enables us to explore self-regulation of 
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thermal pressures. 



4.1. Overall Evolution 

We begin by describing evolution of Model QAl0x2, which has E = 10 Mq pc~^ and p^d = 
0.05 Mq pc~3. Figure [l] displays snapshots for Model QA10x2 at t/torb = 0, 0.1, and 0.2 to 
show early time evolution. The initial gas disk has a Gaussian density profile with scale height 
Hu, = 134 pc and constant temperature, shown in Figure [T]^ a). Since the initial disk is out of 
thermal equilibrium, it rapidly evolves and separates into two phases, with a cold dense layer near 
the disk midplane sandwiched by diffuse warm gas at larger \z\. At the same time, TI develops 
locally, creating numerous cloudlets in the midplane dense layer. The cold midplane slab has a 
surface density of Sc = 7 Mq pc~^ and a typical sound speed Cc = 1 km s~^. The cold slab has 
Toomre stability parameter Qc ~ 0.3 with Jeans length A2d,c = c^/ (GSc) = 33 pc, so that it is quite 
gravitationally unstable. The slab soon fragments gravitationally to form many dense clouds, which 
grow in size and mass by merger with their neighbors. Massive clouds undergo runaway collapse 
as self-gravity dominates the internal pressure, eventually producing stars and SN explosions when 
the density exceeds pcr- The first SN feedback event occurs at about t/torb = 0.1. Figure [T|^6) 
shows formation of dense clouds and the first SN explosion from Model QA10x2. Subsequent SN 
events drive the gas disk into a turbulent state, as seen in Figure [l][^c). 

The kinetic energy associated with expanding shells disperses dense clouds in the midplane, 
and causes the disk to puff up in the vertical direction. Successive stages of gravitational contraction 
and feedback-induced expansion result in quasi-periodic oscillations of the disk thickness. Warm 
gas located ahead of the expanding shells is swept up by shocks and collected into the shells. Pre- 
existing dense gas becomes even denser from shock compression. Ensuing radiative cooling in the 
postshock regions increases the shell density (e.g., Mufson|[l974 McCray et al.|[T975 ). High-density 



expanding shells disintegrate due to a combination of dynamical processes, forming small dense 
cloudlets that subsequently merge together to grow into new dense clouds. These newly formed 
dense clouds collapse internally and create additional stars when their internal density exceeds the 
threshold value, leading to further SN feedback events that repeatedly stir up and restructure the 
surrounding medium. 

Figure [2] plots temporal evolution of the mass fractions of the various gas components, the 
density-weighted vertical scale height 

\ J paxaz J 

and the SFR surface density. The initial changes in the mass fractions and the disk scale height 
shown in Figure [2] reflect early-time thermal response of the gas to the net cooling function. The 
formation of new dense clouds is quickest at the compression phase of the disk oscillation, as 
evidenced by the negative correlation between /gbc and H shown in Figure [2j Within a few tenths 
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of an orbit, the system evolves into a quasi-steady state in the sense that Ssfr, gas fractions, and 
other statistical properties fluctuate but do not systematically change over time. 

Notice that H in Figure [2] shows quasi-periodic oscillations over the entire evolution of Model 
QA10x2, which also produces temporal variations in other physical quantities. The dominant 
timescale is roughly half of the natural vertical oscillation period, ~ 0.5(7r/G/f3sd)^'^'^- The mean 
value and standard deviation of the disk scale height are {H) = 86 pc and /S.H = 12 pc, respectively, 
where the angle brackets () denote a temporal average over 2 < t/torb < 3. When the disk is 
compressed vertically, it produces more dense clouds and hence more active star formation. The 
enhanced radiative and mechanical feedback from star formation then increases the thermal pressure 
and the velocity dispersion of the gas, causing the disk to re-expand. Disk expansion temporarily 
suppresses star formation activity, which then reduces the total pressure and leads to a decrease in 
the disk scale height. At saturation, the mass fraction of the diffuse component in model QA10x2 
has a mean value (/diff) = 0.77 and fluctuation amplitude of A/dig ~ 0.06. The cold, unstable, and 
warm phases amount to fractions (/c) = 0.46, (/„) = 0.22, and {fw) = 0.32, respectively, of the 
diffuse gas mass. The SFR surface density has a mean value (Ssfr) = 1-9 x 10^^ Mq kpc~^ yr~^ 
and standard deviation ASsfr = 4.0 x 10~^ Mq kpc~^ yr~^- Note that ASsfr is small, since in 
evaluating Ssfr we have already time-averaged SN events over tbin = O.Storb (cf. Fig. [2]). 

Figure [3] displays the density structure (including newly formed dense clouds) and velocity field 
around an expanding shell at i/torb = 2.22, well after Model QA10x2 has reached a quasi-steady 
state. The expanding shell, near the center of the simulation box in Figure [3][| a), was created by a 
SN event at t/torb = 2.18. Figure [3]^6), showing a zoomed-in section of the shell, illustrates that 
dense (internal n ~ 10^ - 10^ cm ) clouds form in regions of converging velocity fields, indicated 
as white arrows. The mean velocities of the dense clouds, represented by black arrows, generally 
follow the background converging velocity fields (with an additional random component), suggesting 
that cloud collisions will ensue. The rectangular section marked in Figure ^^h) is enlarged in 
Figure [3]^ c) to show the internal velocity fields of three selected massive dense clouds. The internal 
one-dimensional velocity dispersion in each cloud is ~ 1 km s~^, which is supersonic since the mean 
sound speed inside the dense clouds is ~ 0.5 km s~^. The dense cloud near (x, z) = (—175, —25) pc 
will have a star formation event at a time At = O.Oltorb after this snapshot. 

Figure Q a) shows the distribution of the gas in the n-P plane from Model QAl0x2, averaged 
over t/toj-];, = 2 — 3. The colorbar labels the mass fraction in logarithmic scale. While a large fraction 
of the gas remains close to thermal equilibrium (given by the solid curve) , a non- negligible portion 
is out of thermal equilibrium (~ 18% by mass departs from equilibrium by |AlogP| > 0.15), since 
the gas is continuously disturbed by turbulent motions]^ The thermal equilibrium curve is for the 
time-averaged heating rate; fluctuations AF = 0.16ro relative to the mean value (F) = 0.76Fo 



''The thermal conductivity adopted is somewhat larger than the realistic value, and the numerical diffusion caused 
by large flow speeds also contributes, which may increase the unstable-mass fraction at the expense of the cold gas 
in our models (e.g., Kim et al. 112008 1. 
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displace the equilibrium curve upward and downward. Variations in heating imply that gas can 
be out-of equilibrium with respect to the mean curve (even if instantaneous thermal equilibrium 
holds). Initially after a SN event, some cold gas is converted to the diffuse warm phase, while later 
shock compression and subsequent cooling during later stages of the shell expansion convert some 
warm gas to the cold phase. 

Figures Q 5, c) plot the probability density functions (PDFs) of thermal pressure and number 
density distributions shown in Figure Q a), respectively. Thick and thin lines denote the mass- and 
volume- weighted PDFs. The range of thermal pressure in our models spans more than three orders 
of magnitude, although most of the mass is near the mode of the PDF. The peak value of the 
pressure PDF corresponds to the mean thermal pressure at the midplane. The volume-weighted 
pressure PDF extends toward very small values mainly due to warm gas at high altitude, while 
self-gravitating dense clouds near the midplane occupy the high end of the mass-weighted pressure 
PDF. The mass- weighted density PDF shows the bimodal shape characteristic of the classical two- 
phase ISM (e.g., [Field et al.|1969[[Wolfire et al.|1995HPiontek Ostriker|2004D , although supersonic 
turbulent motions and frequent phase transitions increase the mass fraction in the unstable phase. 



making the peaks less prominent (e.g., 
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Evolution of other models in Series QA is qualitatively similar to that of our standard model. 
One notable trend is that physical quantities exhibit larger-amplitude fluctuations with decreasing 
S. In low-S models where SN events are rare and intermittent, even a single SN explosion stirs 
up the whole simulation domain because there is not enough mass to limit shell expansion. This 
gives rise to large variations in H, which in turn increases the dispersions of Pth and Ssfr, for 
lower-E models. In models with high S, on the other hand, SN events are frequent and spatially 
correlated. Shell expansion is frequently limited by surrounding dense gas and nearby SN shells. 
Consequently, the temporal changes of the disk scale height in these models are less dramatic than 
in low-density models. 

Compared to Series QA, models in Series QB have smaller H, as a result of a more-confining 
vertical gravitational potential (four times larger Psd)- The resulting SFR is correspondingly larger 
in Series QB compared to Series QA. Series S and G also reach quasi-steady states, and their trends 
with increasing/decreasing T, or follow the same patterns as in Series QA and QB. In particular, 
independent increases in either S or (with the other parameter controlled) produce an increase 
in SsFR- The statistical properties of the models vary depending on the input "environmental" 
parameters (i.e. S and Psd), as we shall describe and explore in the remainder of this paper. 
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4.2. Statistical Properties of the Gas 



We have seen in Section [47L] that after a brief transient, our models approach a quasi-steady 
state, which may be thought of as an approximate thermal and dynamical equilibrium (with fluc- 
tuations about the mean). In this subsection, we present the time-averaged values of the physical 
quantities that characterize the thermal and dynamical properties of the gas. These values will be 
used in Section [5] to compare our numerical results with the analytic predictions summarized in 
Section [2l 

In our models, the total pressure at the midplane consists only of the thermal and turbulent 
components since we do not include a magnetic field. We measure these two pressures directly from 
simulation data as 

fz=-Az/2 I P®{n< ?^GBC ) dxdz 
Jz=-Az/2 J @in<nGBC)dxdz 

Jz=-Az/2 J Q(n<nGBc)dxdz 

where @{X) is 1 if the logical argument 'X' is true and otherwise. These definitions give volume- 
weighted averages of pressure for the diffuse component (all gas at n < ncBC = 50 cm~^) at 
the midplane (the horizontal planes z = ±Az/2). Figure [5]^ a) plots as solid and dotted lines the 
midplane thermal and turbulent diffuse-gas pressures, respectively, in Model QA10x2 as functions 
of time. After a quasi- steady state is reached (i/torb > 1)) the mean values are (Pth/^B) = 
1,680 cm-^K and (Pturb/^e) = 5,440 cm'^K, with fluctuation amplitudes APth/ (-Pth) = 0.21 
and APturb/ (Pturb) = 0.52. Since the midplane includes high-velocity injection regions associated 
with SN, there are large spikes in the midplane value of Pturb- The overall fluctuations of Pth and 
Pturb follow the pattern of variations in H due to vertical oscillations, as shown in Figure [2} 

While we can measure midplane pressures in simulations, the most direct observables are 
mass-weighted velocity dispersions. We calculate the mass-weighted vertical turbulent and thermal 
velocity dispersions of the diffuse component using 



/ pv'^Q{n<nQBc)dxdz 
f pQ{n<nQBc)dxdz 



1/2 



1 /2 

J PQ{n<nQBc)dxdz 
f pQ{n<nQBc)dxdz 



(34) 



The rms total velocity dispersion of the diffuse component in the vertical direction is given by 
Cz.diflf = (^zdiff + ^th diff)"'^''^- Figur e [Hjl^ 6 ) displays the time evolution of Vth,d\s and fz,difr in Model 
QA10x2 as solid and dotted lines, respectively. The vertical turbulent velocity dispersion satu- 
rates at (fz,diff) = 6.8 km s~^ with relative fluctuation amplitude Ai^z^diff / (''^z.diff) = 0.30, while 
the thermal component has a smaller mean value (wth,diff) = 3.7 km s~^ and standard deviation 
Ai'th,difr = 0.2 km s~^. Many spikes in fz,difr reflect energy injection events associated with SN 
explosions. Since the shock-heated gas occupies a very small volume only near the midplane, the 
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thermal velocity dispersion fth,diff averaged over the whole domain varies more smoothly than the 
volume-weighted mean Pth averaged only near the midplane. 

Tables [2] and [3] list the mean values and standard deviations of several physical quantities char- 
acterizing the gas disk, for all models. Here and hereafter, we omit angle brackets for convenience; 
all the symbols represent time-averages over i/torb = 2 — 3, unless stated otherwise. Column (1) 
labels each run as in Table [!} In Table [2| Column (2) gives the logarithm of Ssfr in units of 
M0 kpc~^ yr^^- Columns (3) and (4) give the logarithm of -Pth/^B and -Pturb/^B; respectively, in 
units of cm~^ K. Column (5) lists the midplane number density no of hydrogen in units of cm~^ 
defined in analogy with equation (32) but for n rather than P in the integral. Column (6) gives the 
scale height of the diffuse gas -ffdiff = [/ pz^Q{n<nQBc)dxdz/ J p@{n<nGBc)dxdz]^/^ in units of 
pc. 

In Table [sj Columns (2) and (3) give the turbulent and thermal velocity dispersions of the 
diffuse component in units of km s~^, while Column (4) gives the mass- weighted vertical velocity 
dispersion for all the gas = [J {pvl + P)dxdz/ J pdxdz]^^"^ in units of km s""*^. Column (5) lists 
/diff, the fraction of mass in the diffuse component (by definition, all gas at n < ncBC = 50 cm~^ 
is diffuse). In Columns (6) and (7), we hst a = {vl^^^^^ + w^^,diff)/^th,diff and = Vth,diff/4, 
respectively; these parameters are necessary to test the OMLIO theory. Note that (the 
mass fraction of diffuse gas that is warm) since v"^-^ = fwC^ + fcCc and the thermal speed of 
the warm medium is an order of magnitude larger than that of the cold medium Cc- Also note that 
a in Table |3] (based on mass- weighted velocities or pressures averaged over the box) is close but not 
identical to the ratio -Ptot/-fth at the midplane. Finally, Column (8) gives the numerically-measured 
timescale to convert high-density gas into stars, rgF,GBC = (1 ~ /diff)5^/5^SFR in Gyr units; here 
1 — /diff = /gbc is simply defined as the mass fraction at n > uqbc = 50 cm~^. 

Figure [6] plots the mean values of turbulent and total velocity dispersions (a) Vz,diS (b) Cz^diff, 
and (c) (Tz as functions of SgpR for all models except Series R. The mean values over the whole 
set of models shown in Figure [6] are Uz,diff = 6.8 ± 0.6 km s^^, (Jz^diS = 7.7 ± 0.6 km s~^, and 
cJz = 7.0 lb 0.4 km s~^. It is clear that cTz^diff increases slightly as Ssfr increases, while Uz,diff ~ fiz ~ 
7 km s~^ is more-or-less constant in all models (excluding Series R). The slight increase of Uz^diff 
with SsFR is due to an increase of fth,diff with a higher proportion of warm gas at higher Ssfr, 
although thermal speeds (averaged over both warm and cold gas) are lower than turbulent speeds 
for all models except R50 (see Table |3j). 

The nearly constant value of fz,diff) over two orders of magnitude in Ssfr, owes to a balance 
between driving and dissipation for the turbulent momentum (see Section 5.2 for a detailed discus- 
sion). The total vertical velocity dispersion for the whole gaseous medium is also nearly constant 
in all models, ciz ~ 7 km s~^ (see Fig. |6]c). This is because the higher proportion of warm gas 
in the diffuse medium (raising cxz.diff as Ssfr increases) is counterbalanced by a lower proportion 

of the gas being in the diffuse component (which has higher velocity dispersion than the dense, 

1/2 

dynamically- and thermally-cold GBC component) at higher Ssfr. That is, with /^/fj cJz,diff, 
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the larger a^^dis is offset by smaller fdis, for models with higher Ssfr- 

We note that the values of velocity dispersions given in Table [3] and plotted in Figure [6] 
are mass-weighted averages over the entire simulation volume rather than just averages at the 
midplane (which would be Uth,mid = (-Pth/po)^''^ and fz,mid = (-Pturb/po)^''^, where po = lAnipno). 
We report volume- averaged values because these are the closest to direct observables. However, the 
OMLIO theory (and dynamical equilibrium considerations more generally) use midplane values of 
the pressure, which depend on midplane velocity dispersions. We have found Vth,diff /"Wth.mid ~ 1-3 
and Wz.diflf/t^z.mid ~ 1.3 for all models. The reason for this difference is that the gas is somewhat 
differentially stratified, with cold phase preferentially concentrated near the midplane, which makes 
f^th,mid slightly smaller value than Vth.diff- Also, since the gas density and the turbulent dissipation 
rate increase near the midplane, v^, mid is slightly smaller than dis averaged over the whole volume. 

Figure [7] plots the mean values of (a) a, (b) and (c) fwfdis as functions of Ssfr for all 
models except Series R. There is a weak decreasing trend of a with SsfRi but overall a has a small 
range, ~ 3 — 6. The small range of a = {v^^^ + v'^ dis) / '^th dm i™plies that the ratio of turbulent 
to thermal pressure Pturb/-fth = diff/''^th diff diffuse gas is close to constant (for a given 

/rad) over a very large range of Ssfr- The parameter increases as Ssfr increases since a higher 
heating rate increases the warm-gas mass fraction and fw = v^^ diff/^w ~ /to + (1 ~ /w)Cc/c^ ~ fw 
Note that a = 1 + {v^ ^^g/c^)/ fw, so that with ~ Vz,diS ~ 7 km s~^ (see Fig. 6), the decline in 
a ~ 1 + 1/fw from ~ 6 to ~ 3 is just as expected when fw increases from ~ 0.2 to ~ 0.5. The mass 
fraction of warm gas in the whole medium fwfdiS is nearly constant, implying the warmer diffuse 
gas at higher Ssfr is offset by a higher fraction of the medium in a very dense component (here 
defined as n > ncBC = 50 cm~'^). 

For Series R (see Tables [i] and |3| , fth,difr and fw increase as /rad increases (corresponding to 
increasing heating at given Ssfr)- On the other hand, Vz,difT decreases as /^d increases, for the R 
series. Combining these effects, a = 1 + ^^z diff/^th diff decreases by nearly an order of magnitude 
for increasing /^ad in the R series. At large /radi ^^th.diff exceeds fz,difr- On the other hand, (7z,diff 
and fJz decrease only slightly as /rad increases, while Ssfr decreases by a factor ~ 2. Thus, /rad 
appears to affect primarily the energy distribution between thermal and turbulent components that 
results from star formation feedback, together with the proportions of cold and warm gas, for the 
parameter regime we have explored. 



5. Test of the Thermal/Dynamical Equilibrium Model 

5.1. Vertical Dynamical Equilibrium 

Having obtained the statistical properties of the multiphase, turbulent gas from our time- 
dependent numerical simulations, we are now in a position to examine the validity of the assump- 
tions made in OMLIO, and to compare our numerical results with the predictions of the OMLIO 
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analytic theory. 



We first focus on the vertical force balance between (self plus external) gravity and total 
(thermal plus turbulent) pressure. If dynamical equilibrium holds, the total midplane pressure Ptot 
should match the vertical weight of diffuse gas, -Ptot,DE- Taking = l/vr, we rewrite equation ([l]) 
in terms of fdis = Sdis/S and a^^dis = c^ifway^'^ as 



tot,DE 



Jdiff ^ \ /diffj + 



(2 - UsY + 



2 _L ^'^^IdiSPsd 



1.7 X 10^ kB cm-^K /, 



diflf 



10 Mq pc- 



(2 - /diff) + 



(2-/difr)V37(-^ 
V 7 km s~ 



Psd 



0.1 M0 pc- 



1/2' 



(35) 



10 Mq pc- 



-2 



1/2- 



Since the second term in the square brackets of equation (35) dominates for the range of 



parameters we have explored (suitable for outer disks) , an approximate form for equation ( 35 ) is: 
Ptot.DE ~ /diffa,,diffS(2Gpsd)'/' ~ /dVff fT.S(2Gpsd)'/' (36) 



1.0 X 10*A;b cm-^KQ 



1/2 



diff 



7 km s-i/ VlO M0 pc 



Psd 



0.1 Mq pc- 



1/2 



— 1/2 

where we take Cz^diS ^ fdm based on the fact that the velocity dispersions of very dense gas 
are smaller than those of the diffuse component. This is similar to the formula adopted by [Blitz 



Sz Rosolowsky (2004, 2006), except that our expression includes the correction factor fdis and 



allows for the dark matter contribution to psd (see also OMLIO). Although the factor /dig in 



equation (36) is close to unity in outer disks, this correction would be quite important in inner-disk 



regions where gas is dominated by gravitationally-bound GMCs. For the current models, we note 

1/2 

that /diff<Tz,diff/c«; = afdi{ifwCw/(^z,diS ~ /dig'^z/cto ~ 1-0 insensitive to model parameters since 
o'z ~ <7z,diff ~ ~ 7 — 8 km s~^, a ~ 4 — 5, and /diff /to ~ 0.2 — 0.3 if /rad = 1- Thus, if dynamical 
equilibrium is satisfied, we expect the midplane pressure to correlate well with S^/p^. 



Figure |8]^ a) plots the midplane total pressure of the diffuse component Ptot = Pth + -Pturb 
measured from the simulations (as listed in Table [2]) as a function of Tj^p^d for all models. The 
errorbars denote the standard deviations of the pressure fluctuations. The dynamical-equilibrium 



prediction of equation (35) (or the approximation in equation 36) for Ptot.DE can be evaluated 
directly from the model inputs S and psd in Table T] and simulation results for /diff, I'z.diffi and 
t^th.diff listed in Table [sj In the lower panel of Figure [8 a), we plot the relative difference between the 
measured Ptot and Ptot,DE computed from equation ( [35| . These values agree with each other within 
13%. This close agreement verifies that effective hydrostatic equilibrium is indeed satisfied. In 
addition, this suggests that the midplane total pressure in a star-forming disk is set by environmental 
parameters such as the gas surface density, external gravity, the level of the turbulence, etc. Since 
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the parameters appearing in equation (35) can be inferred relatively directly from observables for 
spatially-resolved face-on galaxies (modulo uncertainties in the stellar-disk scale height), the total 
midplane pressure in diffuse gas is an empirically-accessible quantity. 



Adopting the dependence on "environmental" parameters S and psd following equation (36), 
the numerical results are well fitted by 

Ptot = 9.9 X lO^^B cm-^K ( ^ ( ^ . (37) 

VlO Mo pc-2 J VO.I Mo pc-3 J ^ ' 

This fit is overplotted as a dotted line in the upper panel of Figure p[ a). Comparison of the fit to 



the numerical results (equation 37) with the analytic prediction (equation |36[) shows that averaging 

1/2 1 

over our model suite, /diffO'z.diff ~ /difr^z ~ ^-^ km s~ . 

For accounting purposes, we have arbitrarily adopted the choice ncBC = 50 cm~^ as the 
minimum for the dense-gas GBC component. One might be concerned that this may significantly 



affect the value obtained for Ptot.DE- As seen in equation (36), however, Ptot,DE for the present 
models depends on ncBC just through Ptot,DE oc /^j/^ because cTz ~ 7 km s~^ is nearly constant for 
all models. We have checked that if we instead chose ncBC = 100 cm~^, /diff increases by about 
10%, resulting in only about 3% change in Ptot.DE- Thus, for the diffuse-dominated regime studied 
in the present work, i-tot,DE does not depend sensitively on the specific choice for ncBC as long 
as it is large enough. In the regime where gravitationally-bound gas is more important, or where 
self-gravity is comparable to the external gravity, the more exact expression in equation ([T]) (or 



equation 35) should be used for Ptot de- 



While an empirical measure of total midplane pressure can be obtained from spatially-resolved 
observations of S, p^d-, and a^^dm-, pressure-sensitive lines can be used to obtain empirical esti- 
mates of Pth even from unresolved observations. It is thus useful to consider how Pth relates to 
environmental properties in our models. Figure [sjl^ 6) plots the midplane thermal pressure of the 
diffuse component Pth (as listed in Table [2]) as a function of Sy^p^ for all models except Series 
R. The lower panel shows the relative difference between Pth and Pth,DE = ^tot,DE/a as defined 
in equation ([2]), or multiplying equation (35) by 1/a. (Note that this differs slightly from the 



lower panel of Figure |8j^ a) because our measured a is based on volume- averaged rather than mid- 
plane pressures.) The errorbars denote the standard deviations of the pressure fiuctuations. The 
dynamical-equilibrium prediction Pth,DE agrees with the measured Pth at the midplane within 17%, 
excluding Series R. The dotted line in the upper panel of Figure [sF 6) gives our best fit 



Multiplying equation (|36j) by 1/a, the thermal pressure in outer-disk regions is approximately given 
by Pth,DE ~ (/difr/a)o'z,difr5](2G/>sd)^''^- Note that the connection between thermal pressure and 



the parameters S and p^d expressed by equation (38) results from vertical force balance and the 



1/2 

fact that a and /diffO'z^diff ~ /difr'^z nearly constant 
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As seen in Section|4[2| since the amount of energy injected into the thermal component depends 
on /rad, -fth is proportional to fmd for Series R, resulting in significant changes of Pth for the same 
S and psd (see Table [21) . The relation Pth ~ -fth,DE still approximately holds provided that the 



inverse variation of a with /j-ad is included for varying /^ad (see equation 46). Although and 
/diff are insensitive to fj-^d, the large variation of a with /rad implies that the results for Series R 
significantly depart from equation (38). This is why Series R is omitted from Figure [sF 6). 



It is also possible to estimate the midplane density and the scale height. If dynamical equilib- 
rium holds, the mean midplane mass density is given by /Oq.de = Ptot,DE/cr'^^diB = -fth,DE/^^th,diff' °^ 
hydrogen number density no,DE = Po,DE/(l-4mp) by 
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.(39) 



For a Gaussian distribution, the scale height in vertical dynamical equilibrium is 

/diffS 
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(2 - /diff) + 
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1/2- 



(40) 



Figure [9] plots the measured values of (a) the midplane number density uq and (6) the scale 
height of the diffuse gas -f^diff versus the corresponding dynamical-equilibrium estimate given in 
equation (39) and (40), respectively. Our best fits for imposed unity slopes give f^o/"'0,DE = 1-4 
and Hdis/Hdis^BE = 0.87. These differences owe to small differences between the mass- weighted 
thermal velocity dispersion Vth.diff and the slightly- lower midplane value wth,midi as discussed in 
Section [O 



5.2. Thermal Equilibrium and Turbulent Balance 

As described in Section [2j OMLIO hypothesized that the gas disk evolves to a state in which 
both cold and warm phases can coexist at the midplane, at the same thermal pressure, with heating 
balanced by cooling. For given heating rate, a range of pressures between Pmin and Pmax permits 
both a cold and warm phase in thermal equilibrium. For definiteness, OMLIO assumed that the 
midplane thermal pressure Pth in the diffuse medium is comparable to the geometric-mean pressure 

Ptwo — (PminPmax) ^ • 
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In our numerical models, the heating rate evolves with the SFR according to equation (24) 



Assuming the Jpuv.meta contribution is negligible, the geometric-mean pressure is given by equation 



(25), corresponding to Ptwo/^B = 3.1 x 10^ cm ^ K(/i.adSsFR/SsFR,o)) where the coefficient is 
slightly different from that in equation ([s]) since the adopted cooling function in our simulations 



is slightly different from that in Wolfire et al. (2003). For each model, the mean value of Ssfr 



measured from the simulation sets the mean of Ptwo! the mean midplane thermal pressure is also 
measured (see Section 5.1 and Table[2]). Using these measurements. Figure 10 plots Pth/-ftwo as a 



function of Ssfr for all models. The dotted line is our best fit 

= 79 f ^ ^ (41) 

Ptwo VlO-3Mokpc-2yr-V ' ^ ' 

The measured thermal pressure of the diffuse gas is thus smaller than the geometric-mean pressure, 
but only slightly: Pth agrees with Ptwo within ~ 40% for all models, while Pth varies over more 
than two orders of magnitude for our whole suite of models (see Table [2]) . This proves that the 
assumption P^h ~ -ftwo of the OMLIO theory is a reasonable first approximation. 



Using the numerical result given in equation (41), we are now in a position to evaluate the 



thermal yield from feedback r/th defined in equation (11). We find 

/ y \ -0-09 

Our numerical calibration of r/th gives a value ~ 30% lower for the Solar neighborhood than the 
value 1.2/i.ad adopted in OMLIO, and includes a weak decrease of r/th with increasing Ssfr- 

The tendency for r/th to decrease with increasing Ssfr can be understood as follows. Models 
with higher S and SgFR have a larger diffuse-gas density, and hence shorter cooling times, compared 
to models with lower E and Ssfr- In the n-P plane, a shorter cooling time implies that Pth will 
more quickly drop towards Pmin, such that -Pth/-Ptwo will be slightly lower for higher-S, higher-SsFR 
models. Models with lower S have longer cooling times, such that Pth does not drop as quickly 
after heating events, and remains closer to Ptwo- 

Under the assumption that the dynamics of the gas disk has reached a statistical steady state 
(as Figure [5] indicates), the rates of turbulent driving and dissipation must balance each other. For 
mean momentum and mass rri* per supernova, the rate of injection of vertical momentum per unit 
area per unit time to each side of the disk is Pdriv = 0.25(p*/rra*)SsFR, assuming spherical blasts at 
the midplane (OSll). If the injected vertical momentum is preserved until the gas falls back to the 
midplane, the vertical momentum flux across the disk Pturb would be equal to 2Pdriv If) however, 
the injected vertical momentum is dissipated within a vertical crossing time, then Pturb = -Pdriv 
Finally, if the space-time distribution of star formation sites is such that expanding shells collide 
with each other in the vertical direction, then partial cancellation of injected momentum would 

yield Pturb < -Pdriv 



OSll parameterized the uncertainties in dissipation and driving by introducing a factor fp = 
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-PturbZ-Pdriv Here, we use results of our numerical simulations to directly compare the measured tur- 
bulent pressure with the vertical momentum injected by supernovae in our models. We characterize 
the return on mechanical feedback from star formation using the turbulent yield parameter r/turb de- 
fined in equation (12). The parameter fp is related to ?/turb by ?/turb = 3.6/p[(p*/m*)/3000 km s~-^]. 



Figure 11 plots our measurement of the ratio -Pturb /-Pdriv for all models, as a function of Ssfr- 
The dotted line shows our best fit omitting the R series, 

-0.17 



turb 



0.97 



^SFR 



10-3 Mq kpc^^ yr-1 



Our numerical calibration of the mechanical feedback yield is therefore 

-0.17 

??turb = 3.5 



^SFR 



10-3 Mq kpc 



(43) 



(44) 



where we use p^:/m^: = 3000 km s"^ for all models. The numerical result in equation (43) shows 
that /p « 1 provides a good overall estimate; this is also consistent with the results of simulations 
presented in OSll (for the molecule-dominated starburst regime). The numerical result that fp 
(and r/turb) decrease weakly with increasing Ssfr suggests that vertical collisions of shells become 
more important at higher star formation rates, as would be expected. On the other hand, disks 
with lower SsFR suffer somewhat less momentum dissipation because star formation sites are more 
isolated (in space and time), and shells expand into a more rarefied medium. 

As discussed in OSll, the result Pturb ~ -Pdriv is equivalent to having the dissipation time of 



turbulence comparable to the flow crossing time over the largest energy-containing scale (Stone et 



ar||1998 Mac Low et al.|[T998 ), which here is the vertical disk thickness H^is ~ -ffdifr,DE- Feedback 
provides an input momentum per unit time per unit area of ~ -Pdriv ~ ^SFRP*/'m'*- For a dissipation 
time ~ -ffdiff/^Zjdiff) the dissipation rate of vertical momentum in the diffuse ISM, per unit time 



per unit area is ~ Sf^ ^jg/f^diff 
crossing time provided Pf 



Aurb- Thus, driving is balanced by dissipation on a 



turb 



Pdriv ) as in equation (43). 



Combining equations (11) and (12), we have 

Ptot/fes 



103 cm-3 K 



= r] 



SsFR 



10-3 Mq kpc"2 yr-1' 



(45) 



where r] = ?7th + ^turb is the combined yield of thermal and mechanical feedback, with the respective 



contributions given in equations (42) and (44) from our numerical results. Other sources of vertical 



support that are associated with star formation (e.g. radiation pressure, cosmic rays, and magnetic 
fields driven by turbulence) would contribute additional terms to rj. Since rjth and r/turb decrease 
weakly with Ssfr, the increase of Ptot with Ssfr is slightly sublinear. 



Using equations ( 42 ) and ( 44 ) , we obtain an expression for the ratio between total and thermal 



pressure in the diffuse gas: 



a 



1 + 



Vtmh 
Vth 



- 31 - 



This explains the very weak decreasing trend of a with Hsfr for /rad = 1 (see Fig. [7]a). In 
addition, this imphes the value q ~ 5 adopted by OMLIO (based on empirical evidence) is in 
good agreement with the results of numerical simulations (for /rad ~ !)• Similarly, since fw = 
^'tMiff/c' = <diff/(c^«)> ~ [1 + 3.5/-1(Ssfr/10"3 M0 kpc'^ yr-i)-0-08]-i(a,,diff/c»)^ where 
Cz,diff/c«, = 1.1(Ssfr/10~^ M0 kpc~^ yj.-i-)0.04 £qj, ^ models (see Fig.[6]6). This form is consistent 
with the trend for to increase slightly with increasing SgpR, and to increase significantly with 
increasing /rad (see Table [s]). 

Finally, we note that although turbulent energy dominates over thermal energy in equilibrium 
(unless /rad is large), the radiative heating rate exceeds the rate of heating from dissipation of 
turbulent energy, except in far outer disks. The energy input rate ratio is ryth/f?turb times the 
ratio of the turbulent dissipation time (~ -ffdiff/^^z,diff; see Section [6]) to the cooling time (assuming 
thermal equilibrium). In the Solar neighborhood, the cooling time is ~ 1 Myr, whereas the turbulent 
dissipation time is ~ 20 Myr, implying a rate ratio ~ 5. Moving outward in the disk, the radiative- 
to-turbulent heating rate ratio decreases oc ni?diff/^^z,difr) which is oc S for Uz^difi ~ constant. 



6. Star Formation Laws 



In this section, we compare the SFRs obtained in our numerical simulations to SFR formulae 
that are widely used in the literature, both as fitting functions for empirical studies, and as pre- 
scriptions for star formation in numerical models of galaxy formation/evolution. We also introduce 
a new formula that relates Ssfr to the total pressure in the diffuse ISM. This relation follows the 
general form expected when thermal and dynamical equilibrium are both satisfied, and when both 
thermal and turbulent pressure are controlled by feedback from star formation. 



We begin with the orbital time prescription, expressed as Ssfr oc SO (Kennicutt 1998). A 



relationship of this kind is expected if the star formation timescale is proportional to the orbital 
time, which would be true if star formation is governed by large-scale gravitational instabilities 



and the Toomre Q parameter is near its critical value (e.g. Quirk 1972; Wyse Sz Silk 1989 Silk 



1997 Elmegreen 1997 Kim & Ostriker 2001, 2007; McKee & Ostriker 2007). Figure 12 plots the 



mean values of SgpR from our numerical models as a function of SJl. The dotted line is the our 
best fit SsFR = 0.008SO for an imposed unity slope, while the dashed line denotes the empirical 
relation obtained by Kennicutt (1998), Ssfr = 0.017Sr2. The RMS fractional deviation of the 



measurements compared to the fit is 43%. In our simulations, the sites of star formation are 
mainly small-scale dense clouds formed by local thermal and gravitational instabilities, rather than 
very massive clouds formed by large-scale instabilities. Thus, orbital and epicyclic motions do not 
directly control star formation in our models. Rather, the similarity between the behavior of Ssfr 



and Sf] in Figure 12 refiects the correlation of input parameters chosen for our simulations: we set 
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Q CK T, for all models, and since the specific star formation rate increases with S, it also increases 
with fi. 



We next consider Ssfr as a fmiction of S, as shown in Figure [Tsl a). Also plotted as filled 



and empty contours are the recent pixel-by-pixel measurements of Bigiel et al. (2008, 2010) for 



SgpR and S in the regions inside and outside the optical radius, respectively, of nearby spiral and 
dwarf galaxies. Consistent with the observational results for S ^ 10 Mq pc~^, Figure 13 ' a) shows 



that there can be significant variation in EgpR at a given value of S. A single power- law fit to 
the numerical results gives SsFR = 2.2 x Mq kpc"^ yr-i(S/10 Mq pc~2)i-6 (not shown 

in Figure [Tsja), with 33% RMS fractional deviation. Although the power law we find is similar 
to empirical results, our simulations indicate that a single power-law Kennicutt-Schmidt relation 
S oc is not a good fit in outer-galaxy regions where S ^ 10 Mq pc~^ and diffuse atomic gas 
dominates. Close inspection of Figure [l3| a) shows that individually, the QA and QB series each 
follows a relation close to SgpR oc S^, but these relations are vertically offset from each other. The 
reason the QA series has lower SFR than the QB series is that the latter has four times larger psd 
at a given value of S, and the reason both series approximately follow Ssfr oc is that we have 
set psd oc in both series, as we shall discuss below. 

We remark that the current suite of models is not intended to match the full parameter range 
of observed galaxies, but instead to explore the fundamental physical dependence of star formation 
on environmental conditions using carefully controlled numerical models. Nevertheless, Series QA, 
which includes Solar neighborhood conditions and extends to higher and lower S assuming constant 
Q and sq, follows the observed distribution of Ssfr vs. S quite well. At very low gas surface 
density S = 2.5 Mq pc~^, the results from our models have higher Ssfr than much of the observed 
distribution for far outer disks. This is largely because we chose low input values of sq to show the 
effects of stellar gravity clearly in our controlled series of models (lower sq corresponds to higher 
Psd for a given E - see equation 28). Realistic values of sq in far outer disks are likely to be higher 



(see Section 3.3). Higher sq would reduce the vertical gravity and hence reduce Ssfr (following the 
secular trend of decreasing Ssfr with increasing sq = 0.02 to 0.07 to 0.28 from Series S to QB to 
QA at S = 2.5 Mq pc~^). In addition. Series QA, QB, and S fix /rad = Ij whereas /rad is likely to 
increase in far outer disks because of lower shielding where the dust abundance and S are lower (see 



Section 3.2.2). The models of Series R show that Ssfr systematically decreases with increasing 



/rad for fixed S and psd- Thus, the difference between the present model results and observations 
at low S is simply due to differences between model inputs and ambient conditions of gravity and 
shielding in outer galaxies. This emphasizes once again that S alone does not determine SgFR. 

For typical parameters in outer disks, the weight associated with the external (star+dark 
matter) gravity term oc p^^^ exceeds the weight associated with gaseous self-gravity in equation ^ 
(orjsj) for the dynamical-equilibrium diffuse-ISM pressure i-tot,DE) which is equal to the diffuse-ISM 



1/2 A l 

weight. Since the external-gravity dominates, we have -Ptot.DE oc Sp^^ az as in equation (36) (see 

I I I y 1/2 — 

also Figure [8|a), and Ptot.DE oc 7?I1sfr (equation |45j) so that Ssfr oc S/j^^ az/i] for r/ = r/th + Vtmh- 
Since and the yield parameters r/th, Vtmh 8i.re all close to constant (see Fig. [6]and equations 42 
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and 



1 1/2 

44), we expect Ssfr oc Sp^^ . 



1/2 

Figure 13 6) plots results from the simulations for Ssfr vs. Sp^^ , showing a much tighter 



relationship than Ssfr vs. S in Figure [13[ a). Comparing measured values to the fit in equation 



(47) below, the RMS fractional deviation is 24%. This is consistent with recent empirical findings 



(35), adopting CTz = 7 km s 
(equation 42 with /rad = 1 



that star formation is correlated with the stellar, not just the gaseous, content of galactic disks (see 
Section [T]). 

In both panels of Figure 13 we overplot the simultaneous solutions of equations ([5]), (11), and 

a = 5, and isF,GBC = 1-3 Gyr, along with the numerical fit for r/th 
If we instead adopt r/th = 1, the results are quite similar since rj^^ 
is nearly constant. The black dot-dashed curve takes sq = 0.28 as in Series QA, the red dashed 
curve takes sq = 0.07 as in Series QB, and the blue dotted (sq = 0.02) and green long-dashed 
(so = 1.1) curves bracket the overall range of sq for our model suite (see Table [l]). The predicted 
curve for sq = 0.28 (as in Series QA) follows the observations quite well within the optical radius. 
As discussed above, larger values of sq and /rad are likely present in far outer disks, which would 
produce a steeper Ssfr vs. S relation moving to very low S (outside typical optical radii). The 
agreement between numerical models and the simultaneous solution of equations ([5]), (11), and 



(35) confirms the analytic thermal/dynamical equilibrium theory for star formation developed in 



OMLIO. In that work, comparison to individual galaxies shows excellent agreement when both T, 
and in the theory are set from the observations. 

In panel Figure[l3)^&), the black solid line denotes the power-law solution obtained by combining 
equations (11), (38), and (42) (with /^ad = 1) to obtain a prediction for Ssfr: 

1.1 



SsFR = 2.4 X 10- 



Mq kpc~^ yr ^ 



Psd 



0.55 



(47) 



lOMopc-V VO.lM0pc-3^ 
We note that for outer disk regions, the focus of the present models, the approximation f^m ~ 1 



is satisfied, such that the single equation (36) takes the place of the simultaneous solution of 



equations ([s]) and (35). That is, the prediction for outer-disk star formation is independent of 
tsF,GBC- If) rather than using the numerical fit (42) for r]th, we had instead simply adopted a 
constant value of rjth ~ 1, then we would obtain a very similar form to equation (47), except the 



exponent of psd would be 0.5, the exponent of S would be 1, and the coefficient in front would be 

yr- 



2.2 X 10 '^Vth ^0 ^ ^' Small differences between the numerical results and the analytic 



prediction for outer disk regions are due to the fact that some of the idealizations of equation ( 36 ) 
are not fully satisfied in the numerical models. For example, the S series models at S = 15 and 
20 M0 pc~^ have non-negligible gravity from the gas, which increases -Pth.DE above the estimate 



in equation (36), and results in SsFR exceeding the estimate of (47), which neglects the vertical 



gas gravity. Also, we note that the R series, because it has /rad 7^ 1, is not expected to agree with 



equation (47). In fact, members of the R series lie both above and below the prediction, consistent 



with expectations. 



A prescription for star formation commonly used in numerical simulations of galaxy formation 



- 34 - 



and evolution within a cosmological context is to make the star formation timescale proportional 
to the self-gravitation or free-fall time of the gas, cx p~^/^. In the context of disks, it is natural to 
adopt the mean midplane density po ^-s a reference value, so that the SFR surface density would 
be given by 



5] 

SFR = eff(po):i — , (48) 



where tsfi = [37r/(32Gpo)]"^^^ is the free-fall time at the midplane and eff(po) is a star formation 



efficiency per free-fall time at the mean midplane density. Figure 14 'a) plots Ssfr from the numer- 
ical simulations as a function of S/ig^o- The dotted line shows our best fit Ssfr = 0.008(E/tff^o) 
for an imposed unity slope. Note that €s{po) = 0.008 is similar to (but slightly smaller than) the 
value es{pcr) = 0.01 imposed at high density (n cr ~ 500 cm ) for star formation to occur in the 
numerical models. The free-fall time prescription gives a tighter relation than Ssfr vs. Sfi or 
5^SFR vs. S, but there is still scatter. 

Although the free-fall time is commonly adopted as the controlling dynamical timescale, in 
many circumstances self-gravity is less important in confining and condensing gas than the gravity 
of the stars and dark matter. For a given total velocity dispersion (Tz^difTi the vertical dynami- 
cal time is related to the disk thickness by tdyn = ^^diff/o'z,difr- Since -ffdiflf = 5]difr/(\/27r/Jo) = 
Cz,diff/(47rG/9sd)^''^ if external gravity dominates, or H^m = Cz,diff/(''r^G/9o)^^^ if gas self-gravity 
dominates, tdyn ~ 0.3/(Gpmid)^^^ with /Omid = Po + Psd includes both limits. If self-gravity domi- 
nates, tfifi = l.Ttdyii) but if Po <^ Psd, idyn ^ %,0i and the "external" gravity sets tdyn and -f^diff- 

For a disk with significant turbulent contribution to the total velocity dispersion a^^is, ^dyn 
is comparable to the vertical crossing time tver = -f^diflf/'Wz.difr- The vertical crossing time is the 
timescale for turbulence to be dissipated, reducing the disk thickness and raising pQ. For a low 
filling-factor cloudy medium, small, cold, dense clouds can also "fall" to the midplane due to 
the combined vertical gravitational force of stars, dark matter, and gas. When they reach the 
midplane, these small, dense clouds collide dissipatively, collecting into high-mass clouds that are 
internally gravitationally unstable and make stars. The vertical crossing tiniG tver 

is thus expected 

to control how rapidly the diffuse cold component collects into self-gravitating clouds and initiates 
star formation. 



Figure [14[ 6) plots Ssfr from the numerical simulations as a function of S/tver- The dotted 
line indicates our best fit Ssfr = 0.0025(S/tvcr) for an imposed unity slope. The coefficient of this 
fit denotes the star formation efficiency per vertical dynamical time e^cr = 0.0025. The measured 
SFR surface density is well described by the vertical dynamical time prescription, although there 
is still scatter (but slightly less than in Figure [l4)a). The RMS fractional deviations of measured 
SsFR compared to the estimated Ssfr are 26% and 21% for the free-fall time and the vertical 
dynamical time prescriptions, respectively. 



The good correlations shown in Figure 14 for both the tg^ and t^ev prescriptions are presumably 



because both implicitly have similar scaling to Ssfr oc Sy^psd (shown in Figure 13 f*). Since 
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1 /2 

Po ~ Po,DE oc S/ffdiff oc Spg^ (when external gravity dominates), S/tg^o is basically proportional 
to S'^/^Pg^^. For Series QA and QB, S^/^Pg^'^ oc p^^^ because we take psd oc for these models. 
Thus, S/tfr,o oc ^■\/Psd for Series QA and QB. Although Series S and G have somewhat different 
input parameter dependence, the parameter coverage of these model series is not extensive enough 
to reveal a clear difference between Ssfr oc Sy^psj and Ssfr oc E/ifj o- For regions dominated by 
external gravity, we have t^er ~ {^'^Gpsd)~^^'^crz,dis/vz,diS, so that S/t^, 
''^z.difr for our models (and for the real ISM). 



^ ^Pld since (Tz,difr 



We note that the vertical dynamical time prescription for star formation is closely connected 
to the regulation of turbulent pressure by feedback from star formation, and to the relationship be- 
tween input momentum and mean velocity dispersion in the disk (OSll). As shown in Section 5.2 a 



balance between turbulent momentum driving and dissipation is achieved in our models. If Ssfr = 
everS/tver, the momentum driving rate per unit mass becomes 2Pdriv/'^ = 0.5evcrP*/("i*tvcr)- 
Equating this with the expected turbulence dissipation rate ~ 0.5v'^^-^g/Hdis = O.Swz.difr/^ver, we 
obtain fz,diff ~ everP*/m^:. Using our adopted value p^/m^ = 3,000 km s"^ and the efficiency 
ever = 0.0025 measured from our numerical models, this yields fz^diflf = 7.5 km s^^, remarkably 
similar to the mean value fz,difr = 6.8 km s^^ obtained from our numerical simulations. 

As argued in Section [2j energy and momentum feedback from star formation are often the 
dominant sources of heating and turbulence driving, in which case both P^h and Pturb (and therefore 
Ptot) in the diffuse ISM are predicted to vary approximately oc Ssfr (see equations 11 and 12). As 



we show in Section 5.2, our simulations indeed evidence near- linear relations. In Figure 15, we plot 
the measured Ssfr as a function of (a) the measured (Pth/^B)//rad and (c) the measured -Ptot/^B, 
for all of our numerical models. All quantities are time- averaged. Note that the thermal pressure 
is divided by /rad to compensate for the effect of the varying assumed heating efficiency (F/Ssfr). 
The dotted lines in panels (a) and (c) are obtained from equations (11) and (45), respectively, with 



numerical calibrations (42) and (44) for the feedback yields r/th and ?7turb- The dashed line in panel 



(c) plots our best fit omitting the R series: 



SsFR = 2.6 X 10"^ Mq kpc"2 yr^i 



i^tot/A:E 



10^ 



cm" 



-3K 



1.18 



(49) 



The power slightly steeper than unity reflects the weak decline of feedback yields r/th and r/turb with 



5^SFR; as discussed in Section 5.2 (cf. equation 45) . Comparing equation (49) with equation (13), 
we see that our numerical results yield rj = 3.9[Ptot/(10'^fcB cm~^ K)]"''-^^ (for /^ad = 1), quite close 
to the estimate r/ ~ 5 obtained by combining the theory of OMLIO and OSll (see Section |2]). 

In addition to heating/cooling and turbulent driving/dissipation balance, vertical dynamical 
equilibrium is expected to apply, so that the total diffuse-gas pressure at the midplane is equal to 
the vertical weight Ptot = -ftot.DE- Thus, a hallmark of self-regulated star formation, when thermal, 
turbulent, and dynamical equilibrium are all satisfied, is that a relation close to SgpR oc Ptot,DE 
is expected to apply (see equation 13). To the extent that a ~ const., we also expect Ssfr oc 
-fth,DE//rad = ^tot,DE/(o/rad) • In Figure 15 we plot the measured Ssfr from numerical simulations 
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as a function of (6) (i-th,DE/^B)//rad5 and (d) i-tot,DE/^B> for all models. The dynamical-equilibrium 



pressures are computed from input parameters S and psd using equation ( 35 ) and mean measured 
values of /diff, i7z^diff) and a for each modelj^ Dotted and dashed lines are as for Figures 15 ' a,c). 



Figures 15 c,d) show that Ssfr is extremely well correlated with Ptot and -Ptot,DE- The RMS 



fractional deviations of the numerical results from the relation given in equation (49) are only 14% 
and 16% for Ptot and -Ptot,DE, respectively. The correlation is worse if the R series is included. This 
is because 77th oc /j-adj so that higher /i-ad reduces Ssfr compared to other models with the same 
midplane pressure. 

Based on the results of our numerical simulations, we conclude that star formation rates should 
be most closely correlated with the total midplane pressure of the diffuse gas, as in equation 
The relation between diffuse-gas pressure and star formation rate has less scatter than the relation 
between SgpR and the gas surface density S alone, or the combination The relation between 
SgpR and Ptot is also more general than Ssfr oc Sy^/)^ (which applies when external gravity exceeds 
gas self-gravity and (Tz,diff ~ const.), or EgpR oc S/tver (which applies for turbulence-dominated 
disks with Ever ~ const.). In regions dominated by diffuse gas, it is fundamentally the weight of the 
ISM that regulates star formation rates, since star formation rates must adjust until the pressure 
driven by feedback matches this weight. For outer disks that are diffuse-dominated {fdis ~ 1); 
the weight (or Ptot,DE, given by equation [7] or by the approximation in equation [S]) depends only 
on S, Psd) and a^- As noted above, an increase in /j-ad (which would be associated with low dust 
abundance) leads to a decrease in SsFR oc Ptot,DE/??i because r]th oc /rad- 



7. Summary and Discussion 

In this paper, we have used time-dependent numerical simulations to investigate the regulation 
of star formation, as well as the thermal and turbulent properties of the gas, in the regime where 
diffuse atomic gas dominates the multiphase ISM. For the Milky Way and similar galaxies, this 
corresponds to the outer disk - i.e. roughly the Solar circle and beyond. Physical effects included 
in our numerical models (see Section [3]) include differential rotation, Coriolis forces, gaseous self- 
gravity, vertical gravity due to the stellar disk and dark matter halo, interstellar cooling and 



If we compute Ptot.DE from equation (351 using constant values /dift = 0.78 and a^,dis ~ 7.7 km (the mean 
values over the model suite), the best fit to Ssfr vs. Ptot.DB analogous to equation (491 would have a coefficient 
~^ and a power f .05. 



2.2 X fO~ 

^ Although our current numerical models have only explored the diffuse-dominated case, we still expect thermal, 
turbulent, and vertical dynamical equilibrium to hold in the volume-filling diffuse gas even if it is not the dominant 



component of the ISM by mass (see OMLfO, OSll, and Section [2|. In this case, equations (|45| and (|I3l) are still 
expected to hold with near-constant yield coefficients rj, so that Ptot or Ptot.DE would still vary nearly linearly with 
SsFR. It is important to note, however, that in the GBC-dominated case, this is best interpreted as Espr setting /dift 
(by equating 35 and 45 with Espr ~ S/fsF,GBc) rather than the diffuse-ISM weight setting Esfr (see OMLIO). 



If GBCs dominate the mass, Esfr is controlled by the density (and pressure) within the bound clouds. 
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heating, thermal conduction, and feedback from recent star formation in the form of radiative and 
mechanical energy. Although this initial set of models involves a number of simplifications (e.g. we 
consider only a local box representing thin radial- vertical slices so that very large-scale gravitational 
instabilities are absent; we do not include galactic magnetic fields and spiral arms; we omit hot gas 
and treat feedback from SNe via localized momentum injection; we do not explicitly treat radiative 
transfer) , it captures a very important aspect of real ISM disks that is missing in many numerical 
studies of galactic star formation. Namely, the vertical thickness of the disk, and therefore the 
mean gas density, is primarily controlled by (time-variable) turbulence. The turbulent vertical 
velocity dispersion depends on competition between driving by energy inputs from star formation, 
and dissipation through shocks and the mode-coupling turbulent cascade. 

To explore the dependence of Ssfr on environmental parameters, we run models with varying 
total gas surface density S and midplane density p^d of stars plus dark matter. The angular velocity 
Q is set such the Toomre stability parameter Qinit = 2 for a velocity dispersion of 7 km s~^. Our 
models are highly dynamic, but each reaches a statistical steady state within a few tens of Myr. 
In this quasi-steady state, the star formation rate, disk scale height, mass fractions of various gas 
phases, turbulent velocity dispersion, and other physical properties fluctuate about well-defined 
mean values (Fig. [2]). Low-amplitude quasi-periodic oscillations of the disk thickness are correlated 
with episodes of bound cloud formation (at maximum compression) and feedback-driven expansion. 
Small cold clouds repeatedly fall to the midplane and collect (due to self-gravity) into more massive 
clouds, which are then dispersed by feedback from star formation. We use the measured mean 
properties to test the theory of star formation and diffuse-ISM regulation developed in OMLIO and 
OSll, as outlined in Section [2} 

The main results from our simulations are as follows: 

1. We find that most of the gas is at pressures, densities, and temperatures close to thermal 
equilibrium (Fig. [4]). The system evolves to a state in which both warm and cold stable phases 
are present, with mean midplane thermal pressure Pth within ~ 40% of the "two-phase" pressure 



Ptwo = (-Pmin-fmax)^''^, decreasing weakly with increasing SgpR (equation 41). This evolution 
involves continuous re-adjustment of the thermal equilibrium curve, as Ptwo oc F oc Ssfr- Since 
EgpR varies by two orders of magnitude for our model suite, the thermal equilibrium curve shifts up 
and down by the same factor. The midplane thermal pressure increases from Pth/ks ~ 100 cm^^ K 
to ~ 10^ cm^'^K going from low-E, low-psd to high-E, high-psd models (Fig. This finding is 



consistent with the conclusion of Wolfire et al. (2003) that H I should be found in two phases out 



to large distances in the Milky Way (based on an assumed heating rate that declines outward), as 



well as observations indicating both phases are indeed present out to ~ 20 — 25 kpc (Dickey et al. 



2009). Other nearby galaxies also show evidence for both warm and cold atomic gas ( Braun|["l997 



Dickey et al.||20()0 Young et al.||2003 ). The analytic model of OMLIO adopted the assumption that 



the midplane thermal pressure is equal to Ptwo; our numerical results show that this is indeed a 
good first estimate. The result Pth ~ -ftwo implies that radiative heating approximately balances 
cooling. From the point of view of thermal energy replenishment, this means that star formation 
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is highly efficient. 



2. By comparing the prediction of dynamical-equihbrium pressure with the measured time- 
averaged midplane pressure in our numerical simulations, we find that vertical dynamical equilib- 
rium is satisfied within 13% for the total pressure (lower panel of Fig.jsja). For the present models, 
the total weight of the diffuse ISM (-Ptot,DE; given in equation [T| or [35|) is matched by a combination 
of thermal and turbulent pressure. In outer disks, where diffuse gas dominates the total surface 
density Sdifr ~ S, simplified expressions for the total midplane pressure in equilibrium are given by 
equations ([T]) and (|8]). In many outer-disk regions (including the Solar neighborhood), the vertical 
gravity from the stars exceeds that from the gas, such that in equilibrium Ptot oc T,y/psd if the ver- 
tical velocity dispersion is constant. The results from our simulations fit this form well (equation 



37 ) , with a similar result for midplane thermal pressure (equation 38 ) . The numerical results that 
P^-i^ PS y)e and Ptot ~ -ftot,DE demonstrate the validity of the vertical dynamical equilibrium 



assumption in the theory of OMLIO, and confirms prior findings from simulations by Koyama & 



Ostriker (2009b). 



3. Based on our numerical measurements of the thermal and turbulent pressures, we find a ratio 
^tot/^th = a PS 4-5 for essentially all our models (Fig.[7]a) when we fix /rad = (r/ro)(SsFR/SsFR,o)~"^ 
1 (see equation 24). This is consistent with the assumption of OMLIO that a is relatively constant 
for galaxies with shielding properties (and hence Jfuv/^^sfr) similar to the local Milky Way. The 
near-constancy of a results from the fact that both thermal and turbulent pressure are driven by 
feedback (see below). When /rad is varied (for Series R models), corresponding to varying dust 
shielding or FUV heating efficiency, a varies because Pth oc /rad in thermal equilibrium. Higher 
/rad (lower shielding) reduces a following equation (46); for large /rad, -fth can exceed Pturb- 



4. We find that the fraction of diffuse gas in the warm component /„ 



-'th,diff 



/ increases from 



~ 20% to ~ 50% from low- to high-SsFR; when we hold /rad = corresponding to F/Ssfr = const. 
The upper range, with half of the diffuse gas warm (for models similar to the Solar neighborhood). 



is comparable to findings of Heiles & Troland (2003) based on 21 cm emission and absorption 



observations. We find (for Series R) that the warm fraction steeply increases as /rad increases 
(higher F/Ssfr, corresponding to lower shielding by dust). This trend is consistent with the 
finding of Dickey et al. (2000) that the SMC, with a relative metallicity ~ 0.2, has a much higher 



warm-to-cold H I ratio than the Solar neighborhood. We note that in real galaxies, /rad would be 



inversely correlated with S (see Section 3.2.2), which would increase the warm fraction at low T, 
compared to the /rad = 1 models in Series QA, QB, and S presented here. 

5. The time-averaged turbulent vertical velocity dispersions in all of our models are Wz,diff ~ 
7 km s~^, with no systematic dependence on Ssfr (Fig.[6|. Total vertical velocity dispersions az,dis 
in the diffuse medium are larger by ~ 1 — 2 km s~^. The turbulent amplitudes we find, and the 
lack of correlation of Uz^diff with Ssfr, are consistent with observations of H I velocity dispersions 
in the Milky Way and nearby face-on galaxies (Heiles & Troland 2003 Dickey et al. 1990[ [van Zee 



Bryant|[T999{ [Petric Rupen||2007t [Kalberla fc Kerp||2009l ). As discussed in Section [6] (see also 
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OSll), turbulent velocity dispersions v^AiS ~ ^verP*/^* suce expected if the star formation efficiency 
per vertical crossing time is Ever = iver^sFR/S (for tver = -f^diff/i'z.diff), and the momentum injection 
per stellar mass from feedback is p^/m*. Confirming this expectation, the turbulent amplitudes we 
find are consistent with the mean value Ever = 0.0025 measured from our numerical models, for the 
momentum feedback parameter p^,/m^, = 3000 km s~^ used in our simulations. 



6. To assess the balance of turbulent driving and dissipation in our numerical models, we com- 
pare the measured turbulent pressure at the midplane Pturb = Pov^ 
injection rate per unit area Pdriv 



2 (jjg with the fiducial momentum 



0.25(p*/m*)SsFR from star formation feedback. Fig. 11 shows 



that these are approximately equal, decreasing weakly with increasing Ssfr (equation 43). Since 
Pturb represents the characteristic vertical momentum per unit area in the diffuse ISM (Sdiffi^z,diff) 
divided by 2iJcjig/uz dig , this implies the momentum dissipation timescale is comparable to the 
crossing time t^er = -f^diff/^^z.diff) consistent with previous numerical results on turbulent driving 



and dissipation (e.g. Stone et al. 1998 Mac Low et al. 1998). Another way to think of this result 
is that the momentum injected in the diffuse ISM by star formation per unit time is comparable to 
the existing vertical momentum divided by the dynamical time. Thus, from the point of view of 
momentum replenishment, star formation is highly efficient. 

7. We use our numerical models to calibrate the feedback yield parameters ?7th and ?7turb) 
respectively equal to the ratio Pth/^^SFR and Pturb/^^SFR in suitable units (see equations 11 and 



12). Both yield parameters decrease only very weakly with increasing SgpR (see equations 42 and 
44), with thermal yield also depending on the radiation penetration parameter as ryth oc /rad- This 
explains why a = Ptot/Pth = l+'?turb/??th is nearly constant (for /rad = !)• The values ryth ~ 1 x /rad 
and r/turb ~ 4 obtained from our numerical models are consistent with the analytic predictions of 
OMLIO and OSll, respectively. 

8. We compare our numerical results for SgpR to several commonly-used formulae, Ssfr oc Sfi, 

SsFR OC T}' 



SsFR = eff(po)S/%,o (see Figs. 



12 



13 



14). The first two relations are not well 



correlated with the numerical results. The third relation has improved correlation, but this is 
in part because t^^ oc (Gpsd)"^^^ for most of our model suite, and Ssfr is well correlated with 



(Fig. 



13 



see also equation 



47) 



We also compare to the relation Ssfr = CverS/tver for 
tver = Hdis / Vz,diS the Vertical crossing time, which limits how rapidly cold clouds can collect at the 
midplane. The fitted efficiencies are efr(po) = 0.008 and ever = 0.0025, with a stronger correlation 
to the vertical crossing-time than free-fall-time prescription. 

9. The best star formation correlation we find is with the total midplane pressure - either as 
measured in the simulations (Ptot)) or as estimated from vertical dynamical equilibrium (Ptot.DE)- 



Equation (49) fits Ssfr within 16% for all models (excluding Series R), as shown in Fig. 15 Series 



R shows that Ssfr drops if the shielding is reduced (higher /rad)- Our numerical result that Ssfr 
has a near-linear correlation with Ptot.DE is consistent with the analytic models of OMLIO and 
OSll for star formation in diffuse-gas dominated regions - either outer disks or starbursts. The 
near-linear relation between Ssfr and Ptot.DE is also consistent with a similar empirical result 
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Blitz 


2002 


Blitz &: Rosolowsky 


2004 


equation 


45) implies that star format 



found by Leroy et al. (2008), and with the previous empirical findings that molecular gas (the 
immediate precursor of star formation) increases nearly linearly with the ISM pressure ( [Wong fc 

A relationship of the form Ssfr ~ -ftot.DE/"^ (see 
londs to demand: the star formation rate increases 
until the midplane pressure (controlled by thermal and turbulent feedback) balances the vertical 
weight of the diffuse ISM. 

That energy input from massive stars determines the midplane pressure and thus self-regulates 
the star formation rate suggests it is crucial to include stellar feedback, when simulating galactic 



star formation numerically. Indeed, work by Hopkins et al. (2011) contemporary with the present 



study used SPH simulations to show that Ssfr is consistent with the observed Kennicutt-Schmidt 



relations only when feedback is included (see also Dobbs et al. 2011). Without feedback, dense 
clouds collapse in a runaway fashion, increasing the star formation rate by ~ 1 — 2 orders of 



magnitude. Using grid-based simulations of Milky- Way-type galaxies, Tasker ( 2011 ) similarly found 
that star formation rates are at least an order of magnitude higher than observations if feedback is 
not included to drive turbulence and unbind dense clouds that form. Including feedback is known 



to strongly affect the star formation history in long-term simulations of galaxies (e.g. Governato 



et al. 2007). 



Stellar feedback also appears essential for driving and maintaining turbulence in the direction 
perpendicular to the disk plane over many galactic orbits. Other proposed mechanisms for gen- 



erating ISM turbulence include large-scale gravitational instabilities (e.g., Wada et al. 2002 Kim 



et al. 2003 Kim k Ostriker 2007 Agertz et al. 2009 Aumer et al. 2010: Bournaud et al. 2010), 



magnetorotational instabilities (e.g., Kim et al. [2003 Piontek & Ostriker 2004, 2005 2007), and 



non-steady motions generated in spiral shocks (e.g. Kim & Ostriker 2006: Kim et al. 2006, 2010 
Dobbs et al.|2006 > . Turbulence driven by these processes has lower vertical than horizontal velocity 
dispersions, because they all tap galactic rotation. Rotational-gravitational instabilities are able 
to produce turbulence levels comparable to observed values, although vertical dispersions drop to 
^ 4 km s~^ after several galactic orbits (Agertz et al. 2009). In addition, gravitationally-driven 



turbulence is dominated by large scales (i.e., clump-to-clump motions) that do not prevent col- 
lapse within clumps. Without stellar feedback to unbind dense clouds that form, the resulting 
star formation rates are too high. Gravitationally-driven turbulence is likely to be most important 



during the transient, gas-rich early stages of galaxy formation at high redshift (e.g. Ceverino et al 



2010). Characterizing turbulence in galaxies requires subtraction of a "background state," and this 
becomes more difficult to define when there are large secular motions including prominent radial 
flows and collapsing clumps. Even steady spiral shocks create a (steady) azimuthal velocity profile 
that can differ by tens of km s~^ from the background rotation curve. It will be interesting to 
analyze in detail how non-stellar processes combine with stellar feedback to power turbulence over 
both large and small scales in disk galaxies, providing a more complete understanding of galactic 
star formation over all redshifts. 



As noted above, the present numerical models involve radical simplifications compared to the 
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real star-forming ISM. Given the success of these simple models, it is clearly worthwhile to pursue 
further computational modeling along similar lines, improving on the numerical idealizations we 
have made. One of the advantages of local numerical models that resolve ~ pc — kpc is that 
the scales involved directly correspond to those accessible in high-resolution observations of nearby 
galaxies. Results from successive model refinements can be compared to observations to identify a 
"minimal physics" set, incorporating only the most important effects to minimize computational 
cost. 

By employing high-resolution ISM simulations to identify the key processes controlling star 
formation, it will be possible to enhance subgrid models for computational galaxy formation studies 
in the cosmological context. While feedback to drive turbulent pressure plays a dominant role in 
the ISM of the Milky Way and similar galaxies, feedback to drive thermal pressure is likely to 
be increasingly important where there is minimal dust shielding, potentially leading to large /rad 
and ryth ^ ??turb- Some current simulations of dust-poor galaxies at high redshift use subgrid 



shielding models to estimate the abundance of cold, star-forming gas (e.g. Gnedin et al. 2009 



Gnedin &; Kravtsov 2010; Kuhlen et al. 2011). A subgrid model that incorporates both shielding 



and turbulence could potentially bridge over a wide range of redshifts. 



The authors are grateful to the referee for helpful comments on the manuscript. The work 
of C.-G. K. and W.-T. K. was supported by the National Research Foundation of Korea (NRF), 
funded by the Korean government (MEST) under grant No. 2010-0000712. The work of E. C. O. 
was supported by grant AST-0908185 from the U.S. National Science Foundation. 



- 42 - 



REFERENCES 

Agertz, O., Lake, G., Moore, B., et al. 2009, MNRAS, 392, 294 
Agertz, O., Teyssier, R., k Moore, B. 2011, MNRAS, 410, 1391 
Audit, E., & Hennebelle, P. 2005, A&A, 433, 1 
Audit, E., & Hennebelle, P. 2010, A&A, 511, 76 

Aumer, M., Burkert, A., Johansson, P. H., & Genzel, R. 2010, ApJ, 719, 1230 

de Avillez, M. A., & Berry, D. L. 2001, MNRAS, 328, 708 

de Avillez, M. A., k Breitschwerdt, D. 2005, A&A, 436, 585 

Bakes, E. L. O., k Tielens, A. G. G. M. 1994, ApJ, 427, 822 

Basu, S., Mouschovias, T. C., & Paleologou, E. V. 1997, ApJ, 480, L55 

Begelman, M. C., k McKee, C. F. 1990, ApJ, 358, 375 

Bigiel, F., Leroy, A. K., Walter, F., et al. 2008, AJ, 136, 2846 

Bigiel, F., Leroy, A. K., Walter, F., et al. 2010, AJ, 140, 1194 

Bigiel, F., Leroy, A. K., Walter, F., et al. 2011, ApJ, 730, L13 

Blitz, L., k Rosolowsky, E. 2004, ApJ, 612, L29 

Blitz, L., k Rosolowsky, E. 2006, ApJ, 650, 933 

Boissier, S., Prantzos, N., Boselli, A., k Gavazzi, G. 2003, MNRAS, 346, 1215 

Bolatto, A. D., Leroy, A. K., Rosolowsky, E., et al. 2008, ApJ, 686, 948 

Bolatto, A. D., Leroy, A. K., Jameson, K., et al. 2011, ApJ, in press flarXiv:1107.1717p 

Bournaud, F., k Elmegreen, B. G. 2009, ApJ, 694, L158 

Bournaud, F., Elmegreen, B. G., k Elmegreen, D. M. 2007, ApJ, 670, 237 

Bournaud, F., Elmegreen, B. G., Teyssier, R., et al. 2010, MNRAS, 409, 1088 

Braun, R. 1997, ApJ, 484, 637 

Ceverino, D., Dekel, A., k Bournaud, F. 2010, MNRAS, 404, 2151 
Cioffi, D. F., McKee, C. F., k Bertschinger, E. 1988, ApJ, 334, 252 
Daddi, E., Elbaz, D., Walter, F., et al. 2010, ApJ, 714, L118 



-43- 



Dickey, J. M., Hanson, M. M., & Helou, G. 1990, ApJ, 352, 522 
Dickey, J. M. & Lockman, F. J. 1990, ARA&A, 28, 215 

Dickey, J. M., Mebold, U., Stanimirovic, S., & Staveley-Smith, L. 2000, ApJ, 536, 756 

Dickey, J. M., Strasser, S., Gaensler, B. M., et al. 2009, ApJ, 693, 1250 

Dobbs, C. L., & Bonnell, I. A. 2006, MNRAS, 367, 873 

Dobbs, C. L., & Bonnell, I. A. 2008, MNRAS, 385, 1893 

Dobbs, C. L., Bonnell, I. A., Pringle, J. E. 2006, MNRAS, 367, 873 

Dobbs, C. L., Glover, S. C. O., Clark, R C., & Klessen, R. S. 2008, MNRAS, 389, 1097 

Dobbs, C. L., Burkert, A., & Pringle, J. E. 2011, MNRAS, in press f arXiv: 1107.01541 ) 

Elmegreen, B. G. 1997, RMxAC, 6, 165 

Fuchs, B., JahreiB, H., & Flynn, C. 2009, AJ, 137, 266 

Field, G. B. 1965, ApJ, 142, 531 

Field, G. B., Goldsmith, D. W., & Habing, H. J. 1969, ApJ, 155, L149 

Gazol, A., Luis, L., & Kim, J. 2009, ApJ, 693, 656 

Gazol, A., Vazquez-Semadeni, E., & Kim, J. 2005, ApJ, 630, 911 

Genzel, R., Tacconi, L. J., Garcia-Carpio, J., et al. 2010, MNRAS, 407, 2091 

Gnedin, N. Y., & Kravtsov, A. V. 2010, ApJ, 714, 287 

Gnedin, N. Y., Tassis, K., & Kravtsov, A. V. 2009, ApJ, 697, 55 

Governato, F., Willman, B., Mayer, L., et al. 2007, MNRAS, 374, 1479 

Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 

Heiles, C., k Troland, T. H. 2003, ApJ, 586, 1067 

Hennebelle, P. & Audit, E. 2007, A&A, 465, 431 

Heyer, M., Krawczyk, C., Duval, J., & Jackson, J. M. 2009, ApJ, 699, 1092 



Hopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, submitted ( |arXiv:1107M54| ) 

Joung, M. K. R., & Mac Low, M.-M. 2006, ApJ, 653, 1266 

Joung, M. K. R., Mac Low, M.-M., & Bryan, G. L. 2009, ApJ, 704, 137 



- 44 - 



Kalberla, P. M. W. & Kerp, J. 2009, ARA&A, 47, 27 
Katz, N. 1992, ApJ, 391, 502 
Kennicutt, R. C, Jr. 1998, ApJ, 498, 541 

Kim, C.-G., Kim, W.-T., & Ostriker, E. C. 2006, ApJ, 649, L13 

Kim, C.-G., Kim, W.-T., & Ostriker, E. C. 2008, ApJ, 681, 1148 

Kim, C.-G., Kim, W.-T., & Ostriker, E. C. 2010, ApJ, 720, 1454 

Kim, J., Hong, S. S., Ryu, D., & Jones, T. W. 1998, ApJ, 506, L139 

Kim, J., Ryu, D., & Jones, T. W. 2001, ApJ, 557, 464 

Kim, W.-T., & Ostriker, E. C. 2001, ApJ, 559, 70 

Kim, W.-T., & Ostriker, E. C. 2002, ApJ, 570, 132 

Kim, W.-T., & Ostriker, E. C. 2006, ApJ, 646, 213 

Kim, W.-T., & Ostriker, E. C. 2007, ApJ, 660, 1232 

Kim, W.-T., Ostriker, E. C, & Stone, J. M. 2002, ApJ, 581, 1080 

Kim, W.-T., Ostriker, E. C, & Stone, J. M. 2003, ApJ, 599, 1157 

Koo, B.-C, & Kang, J.-H. 2004, MNRAS, 349, 983 

Koyama, H., & Inutsuka, S. 2002, ApJ, 564, L97 

Koyama, H., & Inutsuka, S. 2004, ApJ, 602, L25 

Koyama, H., & Ostriker, E. C. 2009a, ApJ, 693, 1316 

Koyama, H., & Ostriker, E. C. 2009b, ApJ, 693, 1346 

Kroupa, P. 2001, MNRAS, 322, 231 

Krumholz, M. R., McKee, C. F. 2005, ApJ, 630, 250 

Krumholz, M. R., & Tan, J. C. 2007, ApJ, 654, 304 



Kuhlen, M., Krumholz, M., Madau, P., et al. 2011, ApJ, submitted ( |arXiv:1105.2376[ ) 

Lemaster, M. N., & Stone, J. M. 2009, ApJ, 691, 1091 

Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782 

Li, Y., Mac Low, M.-M., & Klessen, R. S. 2005, ApJ, 626, 823 



-45- 



Mac Low, M.-M., & Klessen, R. S. 2004, RvMP, 76, 125 

Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, PhRvL, 80, 2754 
McCray, R., Kafatos, M., & Stein, R. F. 1975, ApJ, 196, 565 
McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565 
Mufson, S. L. 1974, ApJ, 193, 561 

Mouschovias, T. C, Kunz, M. W., k Christie, D. A. 2009, MNRAS, 397, 14 

Ostriker, E. C, McKee, C. F., & Leroy, A. K. 2010, ApJ, 721, 975 (OMLIO) 

Ostriker, E. C, & Shetty, R. 2011, ApJ, 731, 41 (OSll) 

Parravano, A., Hollenbach, D. J., & McKee, C. F. 2003, ApJ, 584, 797 

Parker, E. N. 1953, ApJ, 117, 431 

Petric, A. O., & Rupen, M. P. 2007, AJ, 134, 1952 

Piontek, R. A., & Ostriker, E. C. 2004, ApJ, 601, 905 

Piontek, R. A., & Ostriker, E. C. 2005, ApJ, 629, 849 

Piontek, R. A., & Ostriker, E. C. 2007, ApJ, 663, 183 

Quirk, W. J. 1972, ApJ, 176, L9 

Roman-Duval, J., Jackson, J. M., Heyer, M., et al. 2010, ApJ, 723, 492 

Ryder, S. D., & Dopita, M. A. 1994, ApJ, 430, 142 

Salim, S., Rich, R. M., Chariot, S., et al. 2007, ApJS, 173, 267 

Shetty, R., & Ostriker, E. C. 2006, ApJ, 647, 997 

Shetty, R., & Ostriker, E. C. 2008, ApJ, 684, 978 

Shi, Y., Helou, G., Yan, L., et al. 2011, ApJ, 733, 87 

Silk, J. 1997, ApJ, 481, 703 

Solomon, P. M., Rivolo, A. R., Barrett, J., & Yahil, A. 1987, ApJ, 319, 730 
Sternberg, A., McKcc, C. F., & Wolfire, M. G. 2002, ApJS, 143, 419 
Stone, J. M., k Gardiner, T. A. 2009, NewA, 14, 139 
Stone, J. M., Gardiner, T. A., Teuben, P. et al. 2008, ApJS, 178, 137 



-46- 



Stone, J. M., Ostriker, E. C, & Gammie, C. F. 1998, ApJ, 508, L99 

Tasker, E. J. 2011, ApJ, 730, 11 

Tasker, E. J., & Bryan, G. L. 2006, ApJ, 641, 878 

Tasker, E. J., & Tan, J. C. 2009, ApJ, 700, 358 

Thacker, R. J., & Couchman, H. M. P. 2001, ApJ, 555, L17 

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, 179 

Truelove, J. K., Klein, R. L, McKee, C. F., et al. 1998, ApJ, 495, 821 

van Zee, L., & Bryant, J. 1999, AJ, 118, 2172 

Wada, K. 2008, ApJ, 675, 188 

Wada, K., & Koda, J. 2004, MNRAS, 349, 270 

Wada, K., Meurer, G., & Norman, C. A. 2002, ApJ, 577, 197 

Wada, K., Baba, J., k Saitoh, T. R. 2011, ApJ, 735, 1 

Wada, K., & Norman, C. A. 1999, ApJ, 516, L13 

Wada, K., & Norman, C. A. 2007, ApJ, 660, 276 

Wolfire, M. G., McKee, C. F., HoUenbach, D., et al. 1995, ApJ, 443, 152 

Wolfire, M. G., McKee, C. F., HoUenbach, D., & Tielens, A. G. G. M. 2003, ApJ, 587, 278 

Wong, T., & Blitz, L. 2002, ApJ, 569, 157 

Wyse, R. F. G., k Silk, J. 1989, ApJ, 339, 700 

Young, L. M., van Zee, L., Lo, K. Y., et al. 2003, ApJ, 592, 111 



This preprint was prepared with the AAS IM^jX macros v5.2. 



-47- 



Table 1. Model Parameters 



Model 


E 


Psd 


n 






So 




[ Mo pc-2] 


[ Mo pc-3] 


[ km s"'^ kpc"'^] 


[pc] 


[pc] 




QA02 


2.5 


0.0031 


7 


528 


2048 


0.28 


QA05 


5.0 


0.0125 


14 


269 


1024 


0.28 


QA07 


7.5 


0.0281 


21 


179 


768 


0.28 


QAIO 


10.0 


0.0500 


28 


134 


512 


0.28 


QA15 


15.0 


0.1125 


42 


89 


384 


0.28 


QA20 


20.0 


0.2000 


56 


67 


256 


0.28 


QB02 


2.5 


0.0125 


7 


269 


1024 


0.07 


QB05 


5.0 


0.0500 


14 


134 


768 


0.07 


QB07 


7.5 


0.1125 


21 


89 


512 


0.07 


QBIO 


10.0 


0.2000 


28 


67 


384 


0.07 


QB15 


15.0 


0.4500 


42 


44 


256 


0.07 


S02 


2.5 


0.0500 


7 


134 


768 


0.02 


805 


5.0 


0.0500 


14 


134 


768 


0.07 


S07 


7.5 


0.0500 


21 


134 


512 


0.16 


SIO 


10.0 


0.0500 


28 


134 


512 


0.28 


S15 


15.0 


0.0500 


42 


134 


512 


0.62 


S20 


20.0 


0.0500 


56 


134 


512 


1.10 


G02 


10.0 


0.0250 


28 


190 


768 


0.55 


G05 


10.0 


0.0500 


28 


134 


512 


0.28 


GIO 


10.0 


0.1000 


28 


95 


512 


0.14 


G20 


10.0 


0.2000 


28 


67 


384 


0.07 



Note. — Models S05, SIO, G05, and G20 are identical to QB05, QAIO, QAIO, 
and QBIO models, respectively. All models in Series QA, QB, S, and G have /rad = 
1. Models in the R series (not listed) have the same parameters as model QAIO, 
except /rad = 0.25, 0.5, 2.5, and 5.0 for R02, R05, R25, and R50, respectively. 
All models have = 512 pc except model QA10x2, which is the same as QAIO 
but with Lx = 1024 pc. 
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Table 2. Disk Properties 1 



Model 




lo£ (Ptl, Ik-p.') 


\-* turb/ '^hJ/ 




(Ham) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


QA02 


-4.20 ±0.23 


1.94 ±0.35 


2.61 ± 1.51 


0.05 ±0.08 


342 ± 111 


QA05 


-3.52 ±0.12 


2.53 ±0.22 


3.03 ±0.72 


0.39 ±0.30 


174 ± 37 


QA07 


-3.03 ±0.11 


2.95 ±0.16 


3.57 ± 0.73 


0.82 ±0.46 


132 ± 29 


QAIO 


-2.74 ± 0.11 


3.24 ±0.15 


3.85 ±0.60 


1.12 ±0.58 


92 ±18 


QA10x2 


-2.72 ±0.09 


3.23 ±0.08 


3.74 ± 0.52 


1.32 ±0.40 


94± 10 


QA15 


-2.38 ±0.10 


3.50 ±0.11 


4.08 ±0.51 


1.70 ± 0.69 


70 ±9 


QA20 


-2.06 ±0.10 


3.86 ±0.07 


4.19 ±0.63 


2.76 ±0.70 


51 ± 5 


QB02 


-3.85 ±0.15 


2.29 ±0.23 


2.71 ± 0.63 


0.30 ± 0.26 


157 ± 53 


QB05 


-3.15 ±0.08 


2.81 ± 0.21 


3.45 ±0.61 


0.56 ±0.40 


118 ±38 


QB07 


-2.79 ±0.12 


3.19 ±0.13 


3.75 ±0.71 


1.07 ± 0.55 


77± 20 


QBIO 


-2.58 ±0.09 


3.43 ±0.13 


3.90 ± 0.62 


1.88 ±0.86 


56 ± 12 


QB15 


-2.24 ± 0.05 


3.72 ± 0.09 


4.27 ± 0.55 


2.97 ± 0.98 


44 ±6 


S02 


-3.45 ±0.07 


2.60 ±0.23 


3.25 ±0.56 


0.40 ± 0.30 


110 ±33 


S07 


-2.93 ±0.07 


3.08 ±0.12 


3.62 ±0.68 


1.09 ±0.47 


90 ± 13 


S15 


-2.43 ±0.11 


3.44 ±0.12 


3.96 ±0.70 


1.50 ± 0.94 


95 ±22 


S20 


-2.31 ± 0.10 


3.64 ±0.07 


4.08 ±0.53 


1.79 ±0.58 


87 ± 7 


G02 


-2.82 ±0.08 


3.13 ±0.20 


3.66 ±0.77 


0.88 ± 0.52 


136 ± 34 


GIO 


-2.66 ±0.06 


3.31 ±0.14 


3.86 ±0.60 


1.57 ±0.68 


82 ± 13 


R02 


-2.61 ±0.06 


2.73 ±0.25 


3.97 ± 0.66 


0.98 ±0.54 


99 ± 15 


R05 


-2.72 ±0.08 


2.94 ± 0.15 


3.85 ±0.77 


1.42 ±0.48 


96 ± 17 


R25 


-2.87± 0.12 


3.48 ±0.12 


3.53 ±0.60 


1.11 ± 0.45 


92 ± 12 


R50 


-2.96 ±0.22 


3.69 ±0.16 


3.31 ± 0.37 


1.29 ±0.58 


97± 20 



Note. — The mean values and standard deviations of physical quantities are aver- 
aged over t/tdh = 2 — 3. Col. (2): Logarithmic value of the SFR surface density 
( Mq kpc~^ yr~^). Cols. (3)-(4): Logarithmic values of the midplane thermal and 
turbulent pressures over ( cm^'^K). Col. (5): Midplane number density of hydrogen 
( cm~^). Col. (6): Scale height of the diffuse component ( pc). See Section |4.2| for 
definitions. 
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Table 3. Disk Properties 2 



Model 




V'fth.diff/ 




(/diff) 


a 


fw 


fSF.GBC 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


QA02 


6.20 ± 4.57 


3.27 ± 0.78 


6.86 ±4.67 


0.92 ± 0.07 


4.59 ±5.56 


0.23 ±0.12 


3.31 ± 3.35 


QA05 


6.28 ± 2.95 


3.25 ±0.43 


6.78 ± 2.85 


0.90 ± 0.06 


4.74 ± 3.65 


0.25 ±0.07 


1.73 ± 1.17 


QA07 


6.90 ± 2.78 


3.49 ± 0.46 


7.36 ±2.51 


0.89 ± 0.06 


4.90 ±3.31 


0.29 ±0.08 


0.92 ±0.56 


QAIO 


7.23 ± 2.28 


3.73 ± 0.42 


7.39 ± 2.02 


0.77 ± 0.09 


4.75 ±2.51 


0.32 ±0.07 


1.26 ±0.60 


QA10x2 


6.80 ± 2.07 


3.74 ±0.18 


7.01 ± 1.73 


0.77 ±0.06 


4.31 ± 2.04 


0.32 ±0.03 


1.19 ±0.40 


QA15 


6.95 ± 1.94 


4.05 ±0.30 


7.02 ± 1.71 


0.67 ±0.06 


3.95 ± 1.70 


0.36 ±0.05 


1.20 ± 0.36 


QA20 


6.94 ± 1.67 


4.59 ±0.22 


6.80 ± 1.30 


0.54 ± 0.06 


3.29 ± 1.12 


0.46 ±0.05 


1.06 ±0.27 


QB02 


5.35 ±3.26 


2.93 ±0.39 


5.87 ± 3.18 


0.91 ± 0.08 


4.32 ±4.15 


0.19 ±0.05 


1.64 ± 1.47 


QB05 


7.07 ± 3.34 


3.36 ± 0.45 


7.37 ± 3.02 


0.87 ±0.08 


5.44 ± 4.36 


0.27 ± 0.07 


0.89 ± 0.59 


QB07 


7.15 ± 2.81 


3.62 ±0.31 


7.16 ±2.43 


0.76 ± 0.09 


4.89 ±3.13 


0.30 ±0.05 


1.12 ±0.53 


QBIO 


6.68 ± 2.17 


3.70 ±0.32 


6.55 ± 1.79 


0.67 ±0.08 


4.25 ±2.19 


0.31 ± 0.06 


1.26 ±0.41 


QB15 


7.88 ± 2.05 


4.02 ±0.26 


7.29 ± 1.58 


0.59 ± 0.05 


4.83 ±2.06 


0.35 ±0.04 


1.07± 0.18 


S02 


7.14 ± 3.58 


3.04 ±0.41 


7.25 ± 3.30 


0.84 ± 0.08 


6.53 ±5.75 


0.22 ±0.06 


1.12 ±0.62 


S07 


6.41 ± 2.68 


3.50 ±0.29 


6.80 ±2.44 


0.82 ± 0.05 


4.35 ±2.85 


0.29 ±0.05 


1.12 ±0.38 


S15 


6.76 ± 1.92 


4.25 ±0.31 


7.01 ± 1.76 


0.70 ±0.10 


3.54 ± 1.49 


0.40 ± 0.06 


1.20 ± 0.51 


S20 


6.08 ± 1.43 


4.42 ±0.22 


6.62 ± 1.25 


0.67 ±0.05 


2.89 ±0.91 


0.43 ±0.05 


1.34 ±0.37 


G02 


7.13 ± 2.26 


4.00 ± 0.34 


7.68 ± 2.11 


0.85 ± 0.05 


4.18 ±2.09 


0.37 ± 0.06 


1.01 ±0.39 


GIO 


6.82 ± 1.99 


3.83 ± 0.30 


6.99 ± 1.77 


0.74 ± 0.07 


4.17 ± 1.91 


0.34 ± 0.05 


1.18 ±0.37 


R02 


8.36 ± 2.35 


2.28 ±0.19 


7.85 ± 1.92 


0.78 ± 0.07 


14.51 ± 7.93 


0.11 ±0.02 


0.91 ± 0.32 


R05 


6.80 ± 2.01 


2.87 ± 0.28 


6.84 ± 1.92 


0.80 ± 0.06 


6.62 ±3.50 


0.19 ±0.04 


1.05 ±0.37 


R25 


5.25 ± 1.77 


4.87 ± 0.40 


6.72 ± 1.91 


0.77 ±0.06 


2.16 ±0.80 


0.54 ± 0.09 


1.69 ± 0.65 


R50 


4.49 ± 1.85 


5.70 ± 0.49 


6.78 ± 1.84 


0.85 ± 0.07 


1.62 ±0.52 


0.73 ±0.12 


1.39 ±0.92 



Note. — The mean values and standard deviations of physical quantities are averaged over 4/torb = 2 — 3. 
Cols. (2)-(3): Vertical turbulent and thermal velocity dispersions of the diffuse gas ( km s~^). Col. (4): Total 
vertical velocity dispersion for all gas ( km s~^). Cols. (5)-(7): Mass fraction of the diffuse gas (/diff), the ratio 
of total pressure to turbulent pressure (a), and the square of mass- weighted thermal to wrarm- medium thermal 
speed ('^tij tiiff = /™) ™ the diffuse gas. Col. (8): Timescale to convert dense gas into stars { Gyr). See 
Section 14.21 for definitions. 
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X [pe] 

Fig. 1. — Density snapshots for Model QAl0x2 (logarithmic color scale) at t/toj-\y = 0, 0.1, and 
0.2. The initial single-temperature gas disk (a) evolves rapidly via thermal instability into a con- 
figuration with midplane cold cloudlets sandwiched by outer layers of warm gas. In (6), the first 
SN explosions occur in dense clouds near x = —200 pc produced by mergers and self-gravitating 
contraction of smaller clouds. Subsequent SN explosions disperse the dense clouds and drive the 
disk into a turbulent state (c), in which filamentary structures of cold gas are found at all heights. 
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Fig. 2. - Time evolution in Model QAl0x2 of (a) the mass fractions of the diffuse (/diff, solid) 
and GBC (/gbC) dotted) components, (b) the mass fractions of the cold (/c, dashed), unstable 
(/u, dotted), and warm (/^, solid) phases within the diffuse component, (c) the density- weighted 
vertical scale height H, and (d) the SFR surface density Ssfr- The initial increase of /gbc and 
fc stops at t = O.ltorb = 22 Myr when the first SN event occurs inside a massive dense cloud. 
The model reaches a quasi-steady state after a few tenths of an orbital time, in the sense that the 
physical quantities fluctuate but do not evolve secularly. Note that fdis is positively correlated 
with H. 
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Fig. 3. — (a) Density structure in the whole simulation domain of Model QA10x2 at t/^orb = 2.22, 
including a large, fragmented, expanding shell produced by a recent SN event. (6) The rectangular 
section in (a) is enlarged to identify dense clouds {n > 50 cm~^), outlined by black contours, that 
formed in a region of converging flow where the shell collides with surrounding gas. The white 
arrows represent the background velocity field, while the black arrows show the mean velocity 
of each dense cloud, (c) The section marked in (b) is further enlarged to show internal velocity 
structure of three selected dense clouds. The colorbars (whose range differs from panel to panel) 
indicate number density in logarithmic scale. The sizes of the arrows outside the boxes in (6) and 
(c) correspond to 10 km s~^ and 5 km s~^, respectively. 
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Fig. 4. — (a) Distribution of gas in the n-P/k^ plane for Model QAl0x2, averaged over t/torb = 
2 — 3. The colorbar gives the mass fraction in logarithmic scale. The solid curve marks the locus 
of thermal equilibrium at the mean heating rate of (F) = O.TGFo. Mass- weighted (thick) and 
volume-weighted (thin) probability distribution functions are shown for (b) thermal pressure and 
(c) number density. The vertical dotted lines in (6) and (c) mark the mean midplane thermal 
pressure and number density, respectively, of the diffuse gas. These results show that the system 
evolves to a state in which approximate two-phase thermal equilibrium at a common pressure holds 
for the atomic gas. 
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Fig. 5. — Time evolution of the midplane thermal and turbulent pressures (a) and the thermal and 
turbulent velocity dispersions (b) of the diffuse component for model QA10x2. In (a), P^h initially 

decreases as the gas cools, while Pturb increases rapidly after the gas falls toward midplane and is 
stirred up by SN explosions. After a few cloud formation and feedback cycles (a few 10s of Myr), 
Pth and Pturb reach saturation values of (Pth/^B) ~ 1, 680 cm^^ K and (Pturb/^B) ~ 5, 440 cm~^ K, 
respectively, with a relative fluctuation amplitudes of 0.21 and 0.52. In (6), the velocity dispersions 
saturate at (fth,difr) = 3.7 km and (vz^dis) = 6.8 km s~^, respectively, with relative fluctuation 
amplitudes of 5% and 30%. 
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Fig. 6. — (a) The vertical turbulent velocity dispersion of the diffuse gas Vz^diS, (b) the total 
(turbulent+thermal) velocity dispersion of the diffuse gas (Jz^diS, and (c) the total velocity dispersion 
of all gas (Tz, as functions of the SFR surface density Ssfr for all models except Series R. The points 
and errorbars give the mean and standard deviations over t/torb = 2 — 3. For the whole set of models 
shown in this figure, v^^dis = 6.8ib0.6 km s^^, dz^dis = 7.7ib0.6 km s~^, and = 7.0±0.4 km s^^. 
The dotted lines in all panels show 7 km s^^ for reference. 
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Fig. 7. — Computed values of (a) the ratio of total-to-thermal velocity dispersion for the diffuse 
gas a = 1 + diff/^th diff' (^) square of mass- weighted thermal to warm-medium thermal speed 
/c^ = fw, and (c) the product fwfdis (for fm « /„, the warm gas mass fraction in the diffuse 



v: 



th,difr/''w 

gas and fdis the diffuse mass fraction), as functions of Ssfr for all models except Series R. The 
points and errorbars give the mean and standard deviations over t/torb = 2 — 3. Over more than 
two orders of magnitude in SsfR; the balance between energy input (heating, turbulent driving) 
and energy output (cooling, turbulent dissipation) maintains nearly constant a = Ptot/Ah and 
both warm and cold gas phases. 
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Fig. 8. — Top: Midplane (a) total and (b) thermal pressures of the diffuse gas as functions of S^/Osd- 
The points and errorbars give the mean and standard deviations over t/torb = 2—3. The dotted lines 
in upper panels show fits Ptot/kB = 9-9 x 10^ cm~3j^(E/10 Mq pc-'^){psd/0-l Mq pc-3)V2 and 
Pth/kB = 2.2 X 10^ cm-3K(S/10 M© pc-2)(psd/0.1 Mq pc'^)^/^, respectively. Bottom: Relative 
differences between measured midplane pressures and the dynamical equilibrium estimates using 
equation (35). The mean midplane pressure Ptot varies only 13% relative to Ptot,DE) showing that 
vertical dynamical equilibrium is an excellent approximation. 
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Fig. 9. — Measured versus estimated values of (a) midplane number densities and (6) disk scale 
heights of the diffuse gas. The points and errorbars give the mean and standard deviations over 
*/*orb = 2 — 3. The estimated midplane number density no,DE = /5o,DE/(l-4mp) and scale height 
-f^diff,DE are obtained from dynamical equilibrium as equations (39) and (40), respectively. The 
dashed lines show our best fits no = 1.4no,DE and -ffdiff = 0.87-ffdifr,DE for imposed unity slopes, 
while the dotted lines indicate one-to-one correspondence. 
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Fig. 10. — Measured midplane thermal pressure Pth of the diffuse gas relative to the two-phase 
thermal equilibrium pressure Ptwoj as a function of Ssfr- The points and errorbars give the mean 
and standard deviations over t/torb = 2 — 3 for each model. The dotted line, with a slope of —0.09, 
gives the best fit. Heating/cooling and mass exchange between warm and cold atomic phases enables 
the mean pressure to track the (radiation) energy input from star formation Pth oc Ptwo oc Hsfr 
over more than two orders of magnitude in SsFR- 



- 60 - 



-t3 



O 




T — I — I — r 



T — I — I — r 



I — it — > I— 



-0,5 - 



QA=seri 
QB— seri 
S— seri 
G— seri 
R— seri 



I ' ' ' 

es O 
es A 
es □ 
es ^ 
es m 



0.9?(S,pk/10-Xi^Pc"V"') 



-3,™-lX-0.17 



_I I I L 



-1 1 

log(2srR/lO-'MQkpc-V"') 



Fig. 11. — Measured midplane turbulent pressure -Pturb of the diffuse gas relative to the vertical 
momentum flux injected by star formation Pdriv; as a function of SgpR- The points and errorbars 
give the mean and standard deviations over t/iorb = 2 — 3 for each model. The dotted line with 
a slope of —0.17 gives the best fit. The result Pturb ~ -Pdriv indicates that turbulent driving is 
consistently balanced by dissipation on approximately a vertical crossing time, even though both 
terms vary by more than two orders of magnitude as Ssfr changes. 



- 61 - 



I 



I 

U 
Ph 





o 



-1 F 



-2 - 



3 - 



-4 - 



-5 



I I I I I I I — I — I — I — I I 'I 



■■>>>■■■ 



/; /I ^ 

* T 



QA— series O 
QB— series A 
S— series □ 
G— series 
R-series )(( 



-3 



-2 -1 

log(i:Q) [M^ kpc-' yr-^] 



1 



Fig. 12. — Measured Ssfr as a function of Sil for all models. The points and errorbars give 
mean values and standard deviations over t/torb = 2 — 3. The dotted line shows our best fit 
^SFR = O.OOSSO for an imposed unity slope, and the dashed line shows the empirical result 
SsFR = 0.017Sf^ oflKennicuttlp^). 
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Fig. 13. — SFR surface density Ssfr as a function of (a) S and {b) ^pj^ for all models. The 
points and errorbars give the mean and standard deviations over t/torb = 2 — 3, respectively. In both 
panels, blue dotted, red dashed, black dot-dashed, and green long-dashed lines give the theoretical 
predictions obtained by solving equations ([s]), ( [lT| ), and (35) simultaneously for sq = 0.02, 0.07, 0.28, 
and 1.10, respectively. The parameters cJz = 7 km s~^, a = 5, and tsF.GBC = 1-3 Gyr are held 
fixed for these analytic comparisons, while rjth varies following the numerical fit in equation (42) 
with /rad = 1- Filled and empty contours in (a) show the observational measurements in the 
regions inside ( Bigiel et al.||2008" ) and outside ( Bigiel et al.||2010 ) of the optical radius, respectively, 
for nearby spirals and dwarf galaxies: the contour levels from dark to light correspond to 10%, 
25%, 50%, and 75% of the data. With higher sq and/or /rad at low S (not shown), the models can 
match the observations beyond the optical radius. The black solid line in (b) denotes the power- law 
solution for Ssfr in equation (47). Note that Ssfr is much better correlated with the combination 

1/2 ^ 

Ep^ than with T, alone. 



- 63 - 



-1 



(a) 



& 



O 



-2 - 



-3 - 




QA— series O 
QB-series A 
S— series □ 
G-series 
R-series ¥ 



-3 -1 

log(S/t„.„) [Mq kpc-" yr-'] 



-1 







O 



-3 - 



(b) 




-5 
-2 



QA— series O 
QB— series A 
S— series □ 
G— series 
R-series 5K 



-1 1 

log(£/tver) [Me kpc-^ yr- 



Fig. 14. — Measured SFR surface density T,sfr as a function of (a) S/tff,o and (6) S/tver for all 
numerical models, where tgfi = (37r/(32G/Oo))"^^^ and tver = -ffdifr/^^z,diff are computed using time- 
averaged values of the variables. The points and errorbars give the mean and standard deviations 
over t/torh = 2 — 3. The dotted lines in (a) and (b) show our best fits for imposed unity slopes, 
^SFR = 0.008(S/tff^o) and Ssfr = 0.0025(S/tvei.)i respectively. 
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Fig. 15. — Surface density of star formation Ssfr measured from the simulations as functions 
of (a) the measured midplane thermal pressure of the diffuse gas Pth; (^) the predicted midplane 
thermal pressure Pth.DEi (c) the measured midplane total pressure of the diffuse gas Ptot, and {d) 
the predicted midplane total pressure i-tot,DE- The points and errorbars give the mean and standard 
deviations over t/torb = 2 — 3. Predicted pressures use the dynamical equilibrium equation (35) and 
measured values of /dis, "j and iTz^difi for each model. In (a) and (6) Pth and -Pth,DE are divided by 
/rad to compensate for varying heating efficiency so that Series R may be compared with other series. 
In top and bottom panels, dotted lines are obtained from equations ( 11 ) and (45), respectively, using 
the numerical calibrations ( 42 ) and ( 44 ) . The dashed lines in bottom panels show the best fit given 
by equation (49). The pressures and Ssfr are extremely well correlated, consistent with the idea 
that SsFR adjusts until the pressures (driven by feedback) match equilibrium requirements. 



